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TECHNICAL MEMORANDUM X- 64663 


SUPPLEMENT TO THE ICRPG TURBULENT BOUNDARY 
LAYER NOZZLE ANALYSIS COMPUTER PROGRAM 

INTRODUCTION 


This document serves as a supplement to Reference 1 and contains an 
additional detailed description of the Turbulent Boundary Layer (TBL) Nozzle 
Analysis Computer Program, developed by Pratt & Whitney Aircraft and 
recommended by the Joint Army, Navy, NASA, Air Force (JANNAF) Perform- 
ance Standarization Working Group as a standard reference. 

Among other losses occurring in a rocket thrust ehamber, the viscous 
effects in the boundary layer contribute significantly to the performance 
degradation. It is the objective of this analytical model to determine the _ 
necessary boundary layer parameters which permit the calculation of .the 
thrust deficiency. In order to treat the boundary layer behavior, the edge 
conditions, the wall temperature distribution, and the nozzle geometry must 
be known. Results can only be obtained for a boundary layer development 
along a solid-wall surface area without any film or transpiration coolant flow 
introduction. Heat transfer calculations which are also performed in this 
analysis are not very precise and have only a secondary effect on the perform- 
ance degradation. 

This supplementary document starts with a description of the basic 
calculation sequence. Then the necessary computer program input parameters 
and calculation options are characterized, followed by a summary of the analyt- 
ical results printed by the program. The subsequent section headed by a flow 
chart showing the coordination of all subroutines in the calculation process 
contains a detailed description of each subroutine. After a brief narrative, 
outlining the purpose of the subroutine accompanied by important physical 
equations, the subroutine from a programming aspect is presented. All 
common blocks used in the subroutine, the subroutine calling, the subject 
subroutine, as well as the subroutines called by this subroutine, computer 
library routines, and program built-in library routines, and the calling 
sequence are specified. The underlined parameters in the calling sequence 
represent the terms solved for in the subroutine called. Then the step-by-step 
calculation process in accordance with the computer program listing is reported. 
The symbols used in the. expressions in the subroutine are listed in the front of 
this report. In the appendix a detailed derivation of the important equations 


used in the analysis [1] is presented. This document in connection with 
Reference 1 should provide the user of the subject computer program all 
necessary information to apply the program successfully or to incorporate 
modifications for specific problems. 


ASSUMPTIONS 


The analytical model and associated equations are based upon the 
following assumptions: 


1. The flow is steady, two-dimensional, or axisymmetric. 


2. The boundary layer is confined to a distance from the wall which is 
small compared with either the distance from an axis of symmetry or the 
height of a two-dimensional channel. 

3. The only forces acting on the gas are those caused by pressure 
gradients and skin friction on the wall. 


4. The only changes in total enthalpy are those caused by heat flux 
through the wall. 

5. The flow immediately outside the boundary layer is isentropic and 
one -dimensional parallel to the wall. 1 


6. Pressure is constant through the boundary layer perpendicular to 
the wall. 


7 . The gas follows the perfect gas law and either has constant specific 
heat (calorically perfect) or is thermally perfect [C is a function of temper- 

T P 

ature only and H = f C ( t) dt ] . 


8. The gas has a constant Prandtl number, a viscosity which varies 
as a power of the temperature, and an adiabatic recovery factor equal to 
Prandtl number to the one-third power. 1 


'These assumptions are to be improved. 


2 


9 . The skin friction coefficient is the same as for constant-pressure, 
constant-wall temperature flow on a flat plate at the same free-strcam 
conditions, wall temperature, and momentum thickness. 2 

10 . The Stanton number is the same as for constant-pressure, 
constant- wall temperature flow on a flat plate at the same free-stream. condi- 
tions, wail temperature, and energy and momentum thicknesses. 2 

11 . The Stanton number for unequal momentum and energy thicknesses 

is that for equal thickness multiplied by (<p/ 0) n where n is s small interaction 
exponent. 2 

12. Heat transfer has no effect on skin friction. The term is the 
same as for adiabatic flow. 2 

13. The Stanton number for equal momentum and energy thicknesses 
is related to the skin friction coefficient by von Harman* s form of Reynold' s 
analogy. 


14. Any chemical reactions in the boundary layer affect only the 

driving enthalpy potential for heat flux which can be reflected in the C 
versus T curve. p 

15. The boundary layer shape parameters 0/ 6, A/d, andd*/0 are 
those for l/n-power profiles of velocity and of the difference between stagna- 
tion and wall enthalpy. The exponent n is automatically set equal to seven 
unless otherwise specified in the input. 


CALCULATION PROCESS 


The objective of the analysis is the generation of the necessary param- 
eters to determine the thrust decrement AF_ resulting from the existence 

B » L. 

of a turbulent boundary layer along the nozzle wall. The thrust deficiency 
relationship 


2 These assumptions are to be improved. 
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which is presently not included in the computer program basically represents 
the summation of the shear stress along the nozzle wall. However, this 
summation can be correlated to nozzle exit flow parameters as the derivation 
in the appendix indicates. In the previous equation only the momentum thick- 
ness 0 and displacement thickness 6 * arc representative expressions 
resulting from a boundary layer analysis, while the remaining parameters 
identify potential flow conditions at the boundary layer edge. 


The heat transfer rate to the nozzle wall at any local station, expressed 
by the following relationship 

q = C u p U (H - H ) , 

H e e aw w 


is another important result of this program. Its integrated quantity along the 
wall indicates the energy which is either lost or recovered, depending on the 
cooling concept applied. 


Subsequently, the important equations are listed which must be solved 
within the computer program. 


Momentum equation 
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Energy equation 
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Displacement thickness : 
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Momentum thickness: 
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The power law profiles assuming a 1/7 exponent are used for the 
velocity and enthalpy profiles. 
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for y > 6 


for y ^ A 


for y > A 


With a given nozzle contour r = r(x), the free-stream velocity U 

e 

and density p g are calculated at each axial distance x using the input Mach 
number = M g (x) . The skin friction coefficient and Stanton number C 
are also obtained using empirical relationships. 


Computation of the following parameters 


1. Momentum thickness 

2. Energy thickness 

3. Displacement thickness 

4. Temperature thickness 

5. Velocity thickness 


0 = 0 (x), 

0 - 0 (x) , 

6 * «= 5 * (x) , 
A - A (x), 
5 •= 6 (x), 


6. Temperature or density distribution at the boundary layer edge 

7. Adiabatic wall enthalpy 


H = H (x) 
aw aw 


is performed by using the ideal gas equation 
P = pR T 

and enthalpy relations 


H = / C (t) dt, 

o R 


H = C T, 
P 


for calorically perfect gas, 


or 


V u2 

H = H + (p r ) 7 3 -Ar , 
aw 2 gJ 


— u 

h = h + -At 

o 2 gJ 


H 


/.h- H + (H -H ) fa) -fz 
w o w VA/ 2gJ 


6 


I 
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Sequence of Calculation 

All equations used in the calculation process are presented in Reference 
1. In the appendix of this document a detailed derivation of important equa- 
tions is given, and the equation number is related to the same one used in 
Reference 1. Subsequently, the calculation sequence used in the computer 
program is outlined and schematically shown in Figure 1. 


1. Input values such ns thrust chamber operating conditions and 

stagnation properties (P , T , y , /n , and P ), the nozzle contour r = r (x), 
° oooor 

the free-stream Mach number distribution M ~ M ( x) , and the specific 

G G 


heat versus temperature relationship C 


Cp (t) are necessary require- 


ments. To start the computer program calculation, initial assumptions for 
the momentum and energy thickness must be made, and appropriate option 
indicators must be set. 


2. The following intermediate values M^, 


H e’ 


dH 

-rA U , 

dx e 


dU 

e 

“dx"’ P e’ 


dx ’ P e’ 


dM 
e 

dx ’ 

dH 

w 


dT 

T e’ dx - ’ P e’ C pe’ 


H , H , and — are then internally 
aw w dx 


determined by the program using the associated equations outlined in the 
following section. 

/ A\ 1//? 

3. The important parameter £=(-£> J , varied internally during 

iteration loops, is utilized to calculate the boundary layer thicknesses such as 
A, 6 , 6*, and 6*/0 at a specific location of the nozzle wall. 

4. With empirical relationships the friction coefficient and the 
Stanton number C are determined by iteration. These parameters must 
be known to calculate the heat transfer rate q 


w 


5. Now the Runge -Kutta -Gill method is applied to predict an energy 
and momentum thickness and 0) as a new approximation for a new down- 
stream location x " x + Ax. 




Figure 1* Calculation sequence used in the computer program, 
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(3. The previously outlined computation procedure is repeated using the 
initial approximations from paragraph 5. 

In Figure 1 the previously outlined steps are shown in a diagram format. 
On the left-hand side of each step, the calling subroutine (in brackets), the 
executing subroutine (marked by a solid bubble), and the assisting subroutines 
during the calculation process are outlined. 



Boundary Layer Edge Conditions 


From the chamber stagnation data, the Mach number distribution along 

the nozzle wall and specific heat versus temperature relationship additional 

static flow and thermodynamic properties (y , T , P , U , p J must be gener- 

6 e 6 c 0 

ated to permit the calculation of various integrals as shown in Block 3 of Figure 
1. The local specific heat is determined in subroutine GETPT. At any local 

dM 

0 

station x the Mach number M , its derivative — ; — , and the first derivative of 

e dx 

T , dT /dx are obtained from a quintic spline interpolation routine. In 
e e 

addition, using the definition of enthalpy, 


H = 
o 


H + 
e 



9 


the free-stream velocity is calculated 
U e - H e )3 

where 


T 

e 



Differentiating the previous equation with respect to x, one obtains 

dH dT 

e e 

dx pe dx 


f 

t 
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Differentiating the original enthalpy equation results in 


U 


dU 

e 

c cbc 


dH 

e 

dx 


gJ 


Substituting the right-hand term by the previous expression, the 
following <. quation used in Block 2 of Figure 1 is obtained. 


U 


dV r 
e dx 


dT 

° = - C gJ . 

pc dx 


The derivative of the free-stream density p = p (x) is determined 
by means of the ideal gas equation ~ ° 


P 

e 

P e R T 

e 

If the density change is considered to be very small for a small distance 
Ax along the boundary layer edge, the derivative can be approximated by 

dP e P e (x) - P e (x - Ax) 
dx Ax 

In order to solve the momentum equations in step 5 of the preceding 
section (Sequence of Calculation) , the value of 6 */$ must be given. Equa- 
tions (89) and (119) of Reference 1 are used to determine these quantities. 


Future Improvements 


1. The TBL program uses only the free-stream Mach number M as 

e 

a function of axial nozzle length to describe the boundary layer edge condition. 
All the other important thermodynamic and fluid dynamic properties are calcu- 
lated from the Mach number profile. Furthermore, the molecular weight, 
represented by the specific gas constant, is assumed to remain constant 
throughout the calculation process. 
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It is recommended that the free -stream parameters such as static 
temperature T , static pressure , local density p^ , static enthalpy 

H , local specific heat ratio y y , and the local molecular weight should be 

input to the program which more accurately represents the boundary layer 
edge condition, 

2. The equation representing the thrust decrement caused by the 
boundary layer should be implemented into the program. 

3. The I'elation between the skin friction coefficient and the Reynolds 
number used in the TBL program is based upon low-speed, flow data without a 
pressure gradient. Furthermore, operation of the computer program revealed 
that during numerical iteration the skin friction coefficient showed different 
values when the Reynolds number was greater than 10 4 . In order to avoid this 
problem, it is recommended that the Blasius equation [ 2,3] be used because 

it satisfies experimental low -speed flow data quite well. The Blasius equation 
reads 



0.0128 
“ 05 " 


Re 


6 


with Re 

U 


P U 6 

e_ 

P 


7 


where the density p and the dynamic viscosity p are based upon either the 
free-stream temperature or Eckert' s reference temperature [4,5] . Cebeci 
[6,7] provides another analytical approach to determine the. skin friction 
coefficient, which is compared to some test data. 

For supersonic turbulent flow with a pressure gradient, the following 
modified Blasius equation is recommended 


, _ 0J318 
f “ ReJ 7 ^ ’ 

U 


where 




p m e 

CO CO 


M, 


CO 


This equation describes experimental results obtained by Brott [8] adequately. 


4. The presently used constant power law associated with the velocity 
and enthalpy does not adequately represent the flow condition especially at high 
Mach numbers and for a cooled wall [8-10| . Therefore, it is recommended 
that a variable exponent l/n as a function of the Mach number be incorpo- 
rated into the calculation process, 

5, Presently a constant Prandtl number is input to the program. Since 
the specific heat and the molecular weight change during the expansion process 
in a rocket nozzle, it is recommended that the Prandtl number be internally 
calculated according to the equation 



C 

1 ’ 


C + 
P 





where the specific gas constant 


R 


E. 

31X ’ 


where £ represents the universal gas constant and 3H the mean molecular 
weight. 


TBL COMPUTER PROGRAM DOCUMENTATION 


Tables 1 and 2 describe the TBL computer program input variables 
and the results printed out by the program. Figure 2 shows the TBL computer 
program subroutine linkage. A description of the individual subroutines in the 
computer program is given in the next section. 
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TABLE 1 


TBL PROGRAM INPUT 


Symbol 

Description 

Units 

1. Physical Properties at Stagnation Condition 

P 

o 

PO 

Stagnation pressure 

lbf/ft 2 

T 

o 

TO 

Stagnation temperature 

°R 

y 

o 

GAMO 

Specific heat ratio 
at stagnation condition 


M o 

ZMUO 

Viscosity at stagnation 

lbm/ ( sec -ft) 



condition 


Pr 

PRANDT 

Prandl number 

— 

R 

RBAR 

Gas constant 

( ft-lbf) / ( lbm- 0 R) 


2. 

Contour Geometry r = i 

(x) 

r 

YITAB 

Radius table 

ft 

X 

XITAB 

Axial distance table 

ft 


3. Mach Number in Free Stream M 

= M (x) 




e e 

M 

e 

ZMTAB 

Mach number table 

— 


4. Table of Specific Heat C = 

C (T) 



P 

P 

C 

CPTAB 

Specific heat table at 

Btu/( lbm- 0 R ) 

P 


constant pressure 


T 

TCTAB 

Temperature table 
corresponding to the 
above CPTAB 

°R 

5. Constant Values 

J 

FJ 

Conversion factor 
between thermal and 
work units 

( ft-lbf) /Btu 

g 

G 

Proportionality constant 

( ft/ sec 2 ) ( lbm/lbf) 



inequation F = (M/g)a 
(acceleration of gravity) 
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TABLE 1. TBL PROGRAM INPUT (Continued) 


Symbols 

Description 

Units 

5. Constant Values (Concluded) 


n 

MZ ETA 

Velocity power law jy n 

exponent u/ll = (y/6 ) 

m 

ZMVIS 

Exponent in viscosity- 



temperature relation 

m/p = (t/t ) m 

o o 

n 

ZNSTAN 

Interne .ion exponent in Stanton 
number relation 


6. 

Initial conditions 


e 

xmtial 

d> 

initial 

THETAI 

PHII 

Initial momentum 
thickness 
Initial energy 
thickness 

ft 

ft 

7. 

Wall Temperature (for Option ITWTAB = l) 


T 

w 

TWTAB 

Wall temperature table 

°R 

8. Maximum Length of Step Size 



DXMAX 

Maximum length of 


step size 


9. Tolerances for Iteration Loops 



TOLCFA 

Tolerance in C - C Jt 
iteration loop 

TOLZET 

Tolerance in £ iteration loop 

TOLZMF 

Tolerance used in gas property 
evaluation loops 
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TABLE 1. TBL PROGRAM INPUT (Concluded) 


Option Indicators 


1. Wall Temperature Options 

ITWTAB - -1 T adiabatic wall 

w , . 

temperature (°R) 

0 T ~ constant (°R) 

w 

« 1 T = T (x)(°R) 

w w v ' v ' 


2. Plane or Axisymmetric Flow Option 

EPSZ = 0. Two-dimensional 

= 1. Axisymmetric 


3. Print Option 

IPRINT = 0 Prints only at input values of x 

= 1 Prints at all calculated subintervals 


4. Number of Points in r, x, and M Tables 

1XTAB 4 IXTAB ^ 500 


5. Specific Heat Options 

ICTAB =0 C = constant 

P 

- (3 - ICTAB - 20) Number of points 
in C versus T table 

P 

6. Multiplier to Dimensionalize r and x Tables 

The computer program operates with dimensional quan- 
tities of r and x. If normalized values r and x are 
input, the multiplier SCALE must be input in feet to 
internally calculate dimensionalized values of r and x. 


15 



At every printout point the program prints six groups of quantities 
called for by subroutine DARPRO. 

TABLE 2. TBL PROGRAM OUTPUT 


Symbol 

Description 

Units 

1. Contour Properties 

X 

Axial distance 

ft 

XLARC 

Actual nozzle wall length 
referenced to the first x value 

ft 

YR 

Radius or height of 
contour at X 

ft 

YRP 

Slope of contour 

— 

2. Flow Properties 

ZME 

Mach number 

_ 

TE 

Static temperature 

0 R 

TW 

Wall temperature 

°R 

TAW 

Adiabatic wall temperature 

0 R 

ZMEP 

Mach number gradient 

1/ft 

UE 

Free -stream velocity 

ft/ sec 

PE 

Static pressure 

lbf/ft 2 

3. Boundary Layer 

DELTA 

Velocity thickness 6 

ft 

BDELTA 

Temperature thickness A 

ft 

DELSTR 

Displacement thickness 6* 

ft 

THETA 

Momentum thickness 0 

ft 

PHI 

Energy thickness 0 

ft 

DELSOT 

Shape factor <5*/0 

- 

4. Heat Transfer 

HG 

Heat transfer coefficient h 

S 

Local rate of heat transfer 

Btu/(ft 2 -sec-'R-) 

QW 

Btu/(ft 2 -sec) 

to the wall 

STJMQDA 

Integrated heat transfer Q w 
rate to point X 

Btu/sec 
( axisymmetric) 
Btu/ ( sec-ft) 

( planar) 
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TABLE 2. TBL PROGRAM OUTPUT (Concluded) 


Symbol 

Description 

Units 

4. Heat Transfer (Concluded) 

FORCE 

Drag force in axial F direction 

lbf 

FLAT 

Force normal to X direction F 

y 

for two-dimensional planar flow 
l 

lbf 

5. Internal Integrals 


ZETA, ZI1, 
ZI2, ZI3, 
ZI2P, ZI3P S 
ZI4, ZI5, 
ZI6, ZI7, 
ZI1P 


6. Coefficients 


CF 

Skin friction coefficient C^ 
Stanton number C 

CH 

RTHE 

Reynolds number based on momentum 
thickness R A 

0 

RXLN 

Reynolds number based on wall 
length R^ 

RPHI 

Reynolds number based on energy 

thickness 

<P 

Reynolds number based on 
displacement thickness Rg* 

RDLS 


If the sonic point start procedure is used, the initial values of momen- 
tum thickness, THETAI, and energy thickness, PHII, are printed. 

A new contour table is printed, corrected by displacement thickness 
as follows: 

XCCP Array of corrected normalized axial distance points 

YCCP Array of corrected normal contour points normalized by the 

potential throat radius. 













MAINT B 



Figure 2. TBL computer program subroutine linkage. 
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Description of TBL Program Subroutines 


SUBROUTINE BARCON 

Subroutine BARCON accepts initial conditions and executes a Runge- 
Kutta solution for the boundary layer along the contour. The initial conditions 
are set 

0=0 , 0 = 0 , at x = x . 
o o o 

A step-size Ax, used to increment the axial coordinate x, is determined from 

the input quantity Ax and the local entries in the x-table. Having two 

max ° 

first-order ordinary differential equations 

[lx = f ( X * 0 ’ * ) md S = S ( x * 6 > * ) » 

and using subroutine BARPRO to evaluate the derivatives, a four -term, Runge- 
Kutta numerical solution is used. 

Derivation of corrected axial distance and contour radius is shown in 
the diagram below. 



The corrected contour radius is 
r’ (x) = r(x) - 6* cos /3 , 


I 



where 


c°S ft = gj, 

ds 2 = dr 2 + dx 2 , 





Thus, one obtains 


r f (x) = r (x) 



YCCP . 


The corrected axial distance is 
x’ = x + 6* sin j3 , 

where 


. „ dr dr/ dx 

sm Z 3 - ,j s - (js/dx 


dr/dx 



Thus, one obtains 



COMMON BLOCKS 

COMMON blocks INPUT, INTER, LOOKUP, OUTPUT, and TABLES 
are used. 
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TBL SUBROUTINES 

Subroutine DIRECT calls BARCON. 

BARCON calls subroutines BARPRO, QUITS, START, AND XNTERP. 
FORTRAN SYSTEM ROUTINES 

FORTRAN library routines A LOG and SQRT are used. 

Built-in FORTRAN library routine ABS is used. 

CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL BARCON 
SOLUTION METHOD 

Compute the cubic root of the Prandtl number . 

1 / 

1. PRE 103 = (PRANDT) /3 

Compute term in denominator of the Stanton number. 

2. CHPAR1 = 1.0 - PRANDT + Ln ((6. 0)/(( 5.0) (PRANDT) + 1.0)) 

3. MZETAM = MZETA - 1 

Save the value of exponent in velocity profile. 

4. ZMZETA = MZETA 

5. ZMZETP = ZMZETA +1.0 

6. ZMZETM = ZMETA - 1.0 

7. RMZETA = ( 1. 0)/(ZMZETP) 

8. OOMZET = (1.0)/(ZMZETA) 
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Set the initial value of axial distance, 


9. X = XITAB(l) 

Set the initial value of the difference of axial distances. 
10. DX = 0.0 

Set the initial integrated heat transfer rate. 

12. SUMQDA =0.0 

Set the initial drag force. 

13. FORCE =0.0 

Set the initial normal force. 

14. FLAT = 0. 0 

Set the initial local heat transfer rate. 

15. QW =0.0 

Set the initial heat transfer coefficient. 

16. HG = 0.0 

Set the following parameters. 

17. IXPOS = 1 

18. MX = 0 

19. ITX = 0 

20. IPX = 0 

21. IYX = 0 


22. ITWX = 0 


23. DXIIIIO =■ 0. 0 


24. IBEG - 2 

Initial assumption of adiabatic skin friction coefficient. 

25. CFAGT = 0. 002 

Set the initial value of the sonic point start indicator. 

25. ISTART = 0 

Check whether the sonic point start procedure is to be used. 

27. If THETA I < 0,0, goto 30 
If THETA. I > 0.0, goto 28 

Calculate the shape factor £ based upon initial assumptions. 

28. ZETA = ((PHII)/(THETAI)) RMZETA 

29. Go to 32 

Use the sonic point start procedure. 

30. CALL START 

Set the sonic point start indicator. 

31. ISTART = 1 

32. CFAGP = CFAGT 

Set the energy thickness equal to the input value. 

33. PHI = PHII 

Set the momentum thickness to the input value. 

34. THETA = THETA I 


23 


Set the start value of aNi.nl distance to the Oral value in the axial dis- 
tance table. 


.•if). XIBASK XITAH(l) 

Set the last value of axial distance to the last, value in the. axial distance 
talde. 

3(>. XIKNI) - XITAH(IXTAU) 

Check whether the Mach number table contains at least two values. 

37. If IXTAI5 r; 1 , go to 39 
If IXTAB > 1, go to 38 

38. DXUIIO «= (XITAB( 2) - XIBASE)/( 10. 0) 

Call subroutine BARPRO to obtain C^> C^, dfl/dx, d<f> /dx, q^, h , 
FORCE, SUMQ.DA, XLARC, and so on. 

39. CALL BARPRO(l) 

Call subroutine BARPRO to print the output of BARPRO. 

40. CALL BARPRO(5) 

Call subroutine XNTERP to obtain the radius YR and its first derivative 
YRP, 

41. CALL XNTERP (X, YR, YRP, IYX, XII A B, YITAB, IXTAB, 

CYX, IMX ) 

Save initial displacement thickness and radius of the contour. 

42. DEL - DELSTR 
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45. YMIN c: YR 


44. ONOC - SQIIT (l.O + (YllP) (YRP)) 

Compute the first corrected axial distance point and contour point. 
4 5. XCCP(l) r- X + (DEI jSTR ) ( YR.P) / ( ONOC ) 

4 0. YCCP(l) e YU - ( DEESTlt) /(ONOC) 

Check whether the mach number table contains at least two values. 

47. If IX TAB s? 1, return 
If KTAB > 1, go to 48 

48. Do 101, I = IB EG, IX TAB 
Compute next value of axial distance. 

49. XNEW = XITAB(I) 

50. XMAG = (| XNEW | + |X|)/(2.0) 

51. DXINT = XNEW - X 

52. NX = ( DXINT) /(DXMAX) + 0.99 

53. If NX > 0, go to 55 
If NX < 0, go to 54 

54. NX = 1 

Obtain real value of NX. 

55. ZNX = NX 

Compute "weighted” difference of axial distance values. 

56. DX = (DXINT) /(ZNX) 


25 


57. DXO 2 - (DX)/(2.0) 


58. DXRIIO - (DX)/(10. 0) 

ISO. Do 94, INX el,NX 

Save the value's of energy thickness f and momentum thickness 0. 

(SO. Pi HOLD “ PHI 

01. TIJKOLD = Til 15TA 

Save the value of the axial distance. 

62. XOLD = X 

63. DPIIHtK(l) * (DX) (PH IP) 

64. DTHERK(l) = (DX) (THETAP) 

Compute new value of axial distance. 

65. X = XOLD + DXO 2 

66. Do 80, IRK =2,4 

Check for last time through the Do-loop. 

67. 2! IRK ^ 4, go to 74 
If IRK = 4, go to 68 

Compute new value of axial distance. 

68. X = XOLD + DX 

69. If I (X _ XNEW)/(XMAG)I > 1.0E-6, go to 71 
If | (X - XNE W ) / (XMAG ) I =5 1.0E-6, go to 70 

70. X = KNEW 



Compute new value of energy thickness and momentum thickness. 

71. PHI = PHIOLD + DPHIRK (IRK - l) 

72. THETA = THEOLD + DTHERK (IRK - 1) 

73. Go to 76 

Compute energy thickness and momentum thickness. 

74. PHI = PHIOLD + (DPHIRK(IRK -))(0.50) 

75. THETA = THEOLD + (DTHERK(IRK -l))(0. 50) 

Check whether the energy thickness is negative or zero. 

76. If PHI < 0. 0, go to 85 
If PHI > 0. 0, go to 77 

Check whether the momentum thickness is negative or zero. 

77. If THETA < 0. 0, go to 85 
If THETA > 0.0, go to 78 

Call subroutine BARPRO to obtain d <f> /dx and d0/dx. 

78. CALL BARPRO (IRK) 

79. DPHIRK(IRK) = (DX)(PHIP) 

80. DTHERK(IRK) = (DX) (THETAP) 

Compute energy thickness and momentum thickness according to the 
Runge -Kutta -G ill method . 

81. PHI = PHIOLD + (DPHIRK (l) + (2. 0) (DPHIRK(2) ) 

+ (2.0)(DPHIRK(3)) + DPHIRK (4))/ (6.0) 




82. THETA - THEOLD + (DTHERK(l) + (2. 0) (DTHERK(2) ) 

+ (2.0)(DPHIRK(3)) + DPHERK(4))/(6.0) 

Check whether the energy thickness is negative or zero, 

83. If PHI :s 0. 0, go to 85 
If PHI > 0. 0, go to 84 

Check whether the momentum thickness is negative or zero. 

84. If THETA > 0. 0, go to 88 
If THETA s 0. 0, go to 85 

An error has been made in the calculations; write out an error message 

85. WRITE X, ZME, THETA, PHI 

Call subroutine BARPRO to obtain q^ t h , FORCE, SUMQDA, and 
XLARC. 

86. CALL BARPRO (5) 

87. CALL QUITS 

Call subroutine BARPRO to obtain new C, cL d9/dx, and d^/dx. 

H f 

88. CALL BARPRO (1) 

Select the mini m u m contour radius and its corresponding displacement 
thickness. 

89. IF YR > YMIN, go to 92 
IF YR < YMIN, go to 90 


90. DEL =. DELSTR 



91. YMIN = YR 


Chock for printout of all calculated subintervals. 

92. If IPRINT < 0, go to 94 
If IPRINT > 0, go to 93 

Print data at the calculated subinterval. 

93. CALL BARPRO(5) 

End of INX Do -loop. 

94. Continue 

Call subroutine XNTERP to obtain YR and YRP corresponding to X. 

95. CALL XNTERP (X, YR, YRP, IYX, XITAB, YITAB, EXTAB, 

cyx, imjH - 

96. ONOC = SQRT(l.O + (YRP)(YRP)) 

Compute corrected axial distance points and contour points. 

97. XCCP(I) ^ X + (DELSTR) (YRP) /(ONOC) 

98. YCCP(l) = YR - (DELSTR) /(ONOC) 

Check for printout at input intervals. 

99. If IPRINT > 0, go to 101. 

If IPRINT =£ 0, go to 100, 

Print out data at input intervals. 

100. CALLBARPRO (5) 

End of Mach number table Do -loop. 
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101. Continue 


Compute the potential throat radius. 

102. It POT = YMIN - DEL 

Print out the potential throat radius. 

103. WRITE RPOT 

Normalize the table of corrected contour points, using the potential 
throat radius. 

Nomalized axial distance: 

104. XCCP(l) = (XCCP(l)) /( rpot) 

Nomalized radius: 

105. YCCP(l) = (YCCP(l))/(RPOT) 

106. Do 108, I = IBEG, IXTAB 

107. XCCP(I) = (XCCP(I) )/ (RPOT) 

108. YCCP(I) = (YCCP(I) )./ (RPOT) 

109. WRITE heading for normalized contour point table 
Check whether the sonic point start procedure has been used. 

110. IF 1START ^ 0, go to 113, 

IF ISTART > 0, go to 111 

Print output from sonic point start procedure. 

111. WRITE XCCP(l), YCCP(l) , (l,XCCP(l), YCCP(l), 1= IBEG, 

IXTAB) 


112. Return 


Print regular output. 

113. WRITE (I, XCCP(I), YCCP(I), 1= 1, IXTAB) 

114. Return 

Table 3 gives subroutine BARCON nomenclature. 


TABLE 3. SUBROUTINE BARCON NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

CFAGP 

Guess value of skin friction 
coefficient 


/INTER/, 32 

CFAGT 

Guess value of skin friction 
coefficient 


/INTER/, 25, 
32 

CHPAR1 

Term in the denominator of 
the Stanton number equation 



CYX 

Array of parabola coefficients 
for the nozzle radius table 


/LOOKUP/, 
41, 95 

DEL 

Initial value of displacement 
thickness or that corresponding 
to the minimum nozzle radius 

ft 

42, 90, 102 

DELSTR 

Boundary layer displacement 
thickness 

ft 

/OUTPUT/, 
42, 45, 46, 
90, 97, 98 

DPHIRK 

Array of the first derivative 
of the potential energy thick- 
ness times the axial step size 


DIM, 63, 71, 
74, 79, 81 

DTHERK 

Array of the first derivative 
of the potential momentum 
thickness times the axial step 
size 


DIM, 64, 72, 
75, 80, 82 

DX 

Weighted difference of table 
values of the axial distance 

ft 

/INTER/, 
10, 56-58, 
68, 79, 80 

DXINT 

Difference between axial 
distance and table values of 
axial distance 

ft 

51, 52, 56 

DXMAX 

Maximum length of step size 

ft 

/INPUT/, 52 1 

DX02 

One -ha If DX 

ft 

57, 65 

























TABLE 3. SUBROUTINE BARCON NOMENCLATURE (Continued) 


Symbol 


Description 


Units 


Reference 


DXRHO 


One-tenth the difference 
between axial distances 


ft 


/INTER/ 
23, 38 , 58 


PLAT 


Force normal to x-direc- 
tion for two-dimensional 
planar flow 


lbf 


/OUTPUT/, 

14 


FORCE 


Drag force in axial or 
x-direction 


lbf 

Btu/(ft 2 -sec-°R) 


/OUTPUT/, 

13 


HG 


Heat transfer coefficient 


/OUTPUT/, 

16 


Do-loop counter and 
printout counter 


48, 49, 
106-108, 
111,113 


IBEG 


Sonic line subscript 
counter 


/INTER/, 
24, 48, 
106, 111 


IMX 


Mach number table entry 
indicator and saved 
subscript counter 


/LOOKUP/, 
18, 41, 95 


INX 


Do-loop counter 


59 

/INPUT/, 
92, 99 
/LOOKUP/, 
20 


IPRINT 


Printout indicator 


IPX 


Pressure table entry 
indicator and saved 
subscript counter 


IRK 


Do-loop counter 


66, 67, 71, 
72, 74, 75, 
78-80 


ISTART 


Indicator for sonic point 
start procedure 


26, 31, 110 


ITWX 


Wall temperature table 
entry indicator and saved 
subscript counter 


/LOOKUP/, 

22 


ITX 


1XPOS 


Free-stream temperature 
table entry indicator and 
saved subscript counter 


/LOOKUP/, 

19 


Array position indicator 


/LOOKUP/, 

17 


32 















































TABLE 3. SUBROUTINE BARCON NOMENCLATURE (Continued) 


Symbol 

Description 

IXTAB 

Number of points in the x, 
y, and Mach tables 
(4 IXTAB ^ 500) 

IYX 

Nozzle radius table entry 
indicator and saved 
subscript counter 

mzeta - 

Exponent of velocity 
distribution 

MZETAM 

MZETA minus one 

NX 

Weighting factor for axial 
distance increment 

ONOC 

Square root of one plus the 
slope of the contour squared 

OOMZET 

One divided by Z MZETA 

PHI 

Boundary layer energy 
thickness 

PHII 

Initial value of energy 
thickness 

PHIOLD 

Saved value of energy 
thickness 

PHIP 

Slope of energy thickness 

PRANDT 

Prandtl number 

PRE 103 

Recovery factor - cubic 
root of the Prandtl number 

QW 

Local heat transfer rate to 
the wall 

RMZETA 

One divided by ZMZETP 


Units 


Reference 


/INPUT/, 
36, 37, 41, 
47, 48, 95, 
106, 111, 
113 


/lookup/; 

21, 41, 95 


52-55, 59 


44-46, 

96-98 


/INTER/, 

8 


/OUTPUT/, 
33, 60, 71, 
74, 76, 81, 
83, 85 


/INPUT/, 
28, 33 


60, 71, 74, 
81 


/INTER4 

X 


/OUTPUT/, 

15 











































TABLE 3. SUBROUTINE BARCON NOMENCLATURE (Continued) 


Symbol 


RPOT 


SUMQI3A 


THEOLD 


THETA 


THETAI 


THETAP 


X 


XCCP 


XIBASE 


XIEND 


XITAB 


XLARC 


XMAG 


Description 


Potential throat radius 
corrected for displacement 
thickness 


Integrated heat transfer 
rate to point X 


Saved value of momentum 
thickness 


Boundary layer momentum 
thickness 


Input value of momentum 
thickness; if = - 1, 
sonic point start procedure 
will be used 


Slope of momentum 
thickness 


Axial distance in 
monotonically increasing 
order 


Array of normalized 
corrected axial distance 
points 


First value in axial 
distance table 


Last value in axial 
distance table 


Table of EXTAB values 
(axial distance, x) in 
monotonically increasing 
order 


Arc length of the contour 
up to point x 


One -ha If the sum of 
XNEW and X 


Units 


ft 


Btu/sec( axial) 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


ft 


Reference 


102-105, 
107, 108 


/OUTPUT/, 

12 


61, 72, 75, 
82 


/OUTPUT/, 
34, 61, 72, 
75, 77, 82, 
84, 85 


/INPUT/, 
27, 28, 34 


/INTER/, 
64, 80 


/OUTPUT/, 
9, 41, 45, 
50, 51, 62, 
65, 68-70, 
85, 95, 97 


DIM, 45, 
97, 104, 
107, 111, 
113 


/INTER/, 
35, 38 


/INTER/, 

36 


/TABLES/, 
9, 35, 36, 
38, 41, 49, 
63, 64, 95 


/OUTPUT/, 

11 

50, 69 
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TABLE 3. SUBROUTINE BARCON NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

XNEW 

Successive values of axial 
distance from the XITAB 
table 

ft 

49-51, 69, 70 

XOLD 

Saved value of axial 
distance 

ft 

62, 65, 68 

YCCP 

Array of corrected 
contour points normalized 
by the potential throat 
radius 


DIM, 46, 98, 
105, 108, 
111, 113 

YITAB 

Nozzle contour radius 
table related to IXTAB 
array 

ft 

/TABLES/, 
41, 95 

YMIN 

Saved value of initial or 
minimum radius or height 
of contour 

ft 

43, 89, 91, 
102 

YR 

Nozzle radius or contour 
height 

ft 

/OUTPUT/, 
41, 43, 46, 
89, 91, 95, 
98 

YRP 

Slope of contour 


41, 44, 45, 
95-97 

ZETA 

Shape factor 


/OUTPUT/, 

28 

ZME 

Mach number 


/OUTPUT/, 

85 

ZMZETA 

Real value of MZETA 


/INTER/, 
4-6. 8 

ZMZETM 

ZMZETA minus one 


/INTER/, 6 

ZMZETP 

ZMZETA plus one 


/INTER/, 5, 
7 

ZNX 

Real value of NX 


55,. 56 


i 
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SUBROUTINE BARPRO 


Subroutine BARPRO calculates the boundary layer parameters from 
initial or previously determined energy thickness <p f momentum thickness 0 , 
dfl d0 

and and -j£- along the contour at a point x. 


The interpolation routine XNT&RP is used to evaluate inviscid flow and 

dM a dT 

contour properties and derivatives at point x. such as M , — £. t . r 

d r e’ dx e’ dx’ ’ 

and dx ‘ Then subroutine ZETA.IT is called and the shape factor £ and bound- 
ary layer thicknesses A, 6 , and 6 * are computed. An iteration procedure is 
used to calculate the skin friction coefficient C f and the Stanton number C . 
This procedure is as follows: 1 H 


1. An initial guess C. is made. 

fag 


2 . 


3. 


4. 


The term 




T \ I 
aw 


Subroutine CFEVAL is used 


The terms 


( T Iw) 


and C. = 
fa 


-m 

CfagRg is calculated, 
to evaluate C { asa function of C 

) a re .calculated. 
' 



5. A relative error comparison is made. If 


C 


fa ^fag 

C, 

fag 


TOLCFA, 


convergence is satisfactory. Otherwise a form of Wegstein • s method is used 
to calculate a new guess C^; and steps 2 through 5 are repeated up to a 

maximum of 50 iterations. 
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The derivatives ~ and t- urn now evaluated, from the differential 
clx dx 

equations, at the point x. For the point x at the end of a Runge-Kutta step 
(IND " l) , total heat transfer and drag' quantitic s are also computed. 

This subroutine also provides the data for the printout if called for by 
the input parameter IPRINT. 

The values obtained in BARPRO are as follows: 


dM < (IT (Ui (i 

M e’ 'lit ’ T e’ lit’ P e ! C po» H o* “ ’ Y e ’ 


dU dp d(p U-) 

C o TT ^ ^ . . T I rr\ 

V "dx ’ P e’ dx ’ P e o’ J " ’ P "’ 


dx 


e aw aw 


II* ^2» ^3» *4> ^5* ^G» ^7* ^*1» ^2> ^ 3» 

* c t* — “ d0 d<j> dr 

t,, A, 5, 6 , C f , C H , dx , dx * 

V V V R e’ R x’ R U V R 5*’ F x’ F y* 


COMMON BLOCKS 


COMMON blocks COFIIF, CSEVAL, INPUT, INTER, LOOKUP, 
OUTPUT, SAVED, and TABLES are used. 


TBL SUBROUTINES 

Subroutine BARCON calls BARPRO. 

BARPRO calls subroutines CFEVAL, CETPT, QUITS, SEVAL, 
XNTERP, and ZETAIT. 

FORTRAN SYSTEM ROUTINES 

FORTRAN library function SQRT is used. 

Built-in FORTRAN function ABS is used. 
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CALLING SEQUENCE 

The: subroutine calling sequence is.> 

CALLBARPRO (IND) 
where 

IND - program loop control indicator. 

SOLUTION METHOD 

Determine which program option is desired. 

1. If IND = 1, go to 2 
If IND = 2, go to 2 
If IND = 3, go to 58 
If IND = 4, go to 2 
If IND = 5, go to 150 

Call subroutine XNTERP to obtain the free-stream Mach number M 
and its gradient dM^/dx for a given x. 1 

2. CALL XNTERP (X, ZME, ZMEP , IMX, XITAB, ZMTAB, 

IXTAB, CMX, IXPOS) 

Save IMX. 

3. IXPOS = IMX 

Call subroutine XNTERP to obtain the free-stream temperature 

T and its gradient dT /dx. 
e e 

4. CALL-XNTERP (X, TE, TEP, ITX, XITAB, TITAB, IX TAB, 

CTX, IXPOS) 
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Call subroutine GETPT to obtain the free-stream pressure by 

using M and T . 
e e 

5. CALL GETPT (ZME, PE, TE) 

Call subroutine SEVAL to obtain the specific heat C and enthalpy 
H g by using T^. 

6. CALL SEVAL (l, TE, CPE, HE) 

/ dH e dT e \ 

Determine the enthalpy gradient I = C pe “djT * gJ ) * 

7. HEP = (FJG) (CPE) (TEP) 


Compute the specific heat ratio 



8. GAME = (CPE)/(CPE - ROJ) 


Obtain the dynamic enthalpy (UV2) . 


9. UE202 = HO - HE 


The free-stream velocity U is obtained from equation (9) . 

e 

10. UE = SQRT ((2.)(UE202)) 


Determine the velocity gradient 




11. UEP = -(HEP)/(UE) 


Obtain the free-stream density p 

e 




12. RHOE = (PE)/(TE)/(RBAR) 
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13. If DXRHO * 0.0, go to 10 
If DXRHO = 0.0, go to 14 

Save the initial density gradient. 

14. RHOEP =0.0 

15. Go to 32 

16. If X > XIBASE, go to 20 
If X ^ XIBASE, go to 17 

Save the density p^. 

17. Z1 = RHOE 
Put one as: 

18. Z3 = 1.0 

19. Go to 24 

Call subroutine XNTERP to determine M g = ( ERASE l) 
and dM g /dx = (ERASE2) for ( X -DXRHO) . 

20. CALL XNTERP (X-DXRHO, ERASE1, ERASE2, IMX, XITAB, 

ZMTAB, IXTAB, CMX, IXPOS) 

Call subroutine GETPT to obtain P = (Z4) and T = (Z5) 

G 6 

corresponding to = (ERASEl). 

21. CALL GETPT (ERASEl, Z4, Z5) 

/ P 

Compute the density p 4 corresponding to M . ( p = — - 

G GIG K 1 

\ G 


22. Z1 = (Z4)/(Z5)/(RBAR) 



t 


L- 




Set point one as: 

23. Z3 = 0.5 

24. If X < XIEND, go to 28 
If X ^ XIEND, go to 25 

Set 

25. Z2 = RHOE 

26. Z3 = 1.0 

27. Go to 31 

Call subroutine XNTERP to obtain M = (ERASEl) and 

e 

dM /dx = (ERASE2) corresponding to X + DXRHO. 
e 

28. CALL XNTERP (X + DXRHO, ERASEl, ERASE2, IMX, XITAB, 

ZMTAB, IXTAB, CMX, EXPOS) 

Cal) subroutine GETPT to obtain P = (Z4) and T = (Z5) 

e e 

using the above M g = (ERASEl). 

29.. CALL GETPT (ERASEl, Z4, Z5) 

Compute the density p = (Z2). (p = F /RT ) 

6 6 6 6 

30. Z2-F. (Z4)/(Z5)/(RBAR) 

Approximate the density gradient. 

31. RHOEP = ( (Z2 - Zl)/(DXRHO))(Z3) 

Obtain the mass flow density p U . 

e e 

32. RHOUE = (RHOE) (UE) 
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ft' 


The first derivative of mass flow density: 


iiP c V J 

dx 


33. RHOUEP = (RHOE) (UEP) + (UE) (RHOEP) 


/ T e\ m 

Evaluate the viscosity = M I I 


34. ZMU = (ZMUO) ((TE)/(T0)) ZMVIS 

1/ u 2 

Compute the adiabatic wall enthalpy H = H + Pr % — 

aw e 2 

35. HAW = (HE) + (PRE103)(UE202) 

Call subroutine SEVAL to determine T and C 

aw pw 

= ( ERASE3) using the known H 

aw 

36. CALL SEVAL (2, TAW , ERASE3 , HAW) 

ITWTAB = - 1: adiabatic wall temperature 

= 0: constant wall temperature 

- 1: input wall temperature (variable) 

37. If ITWTAB < 0, go to 38 
If ITWTAB = 0, go to 42 
If ITWTAB > 0, go to 46 

Consider the case of adiabatic wall temperature. 

Set 


38. TW = TAW 

39. HW = HAW 
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The enthalpy gradient along the walls 

40. IIWP = HEP + (PRE103) (UE) (UEP) 

41. Go to 49 

Consider the case of constant wall temperature. 

42. TW = TWTAB(l) 

The wall enthalpy is calculated at 86 of BARSET; 

43. HW = TWTAB (2) 

The gradient of wall enthalpy is zero for constant wall 
temperature option. 

44. HWP =0.0 

45. Go to 49 

Call subroutine XNTERP to find T and obtain dT /dx 

w w 

in the case of variable wall temperature option. 

46. CALL XNTERP (X, TW, TWP , ITWX, XITAB, TWTAB, IXTAB, 

CTWX, IXPOS) 

Call subroutine SEVAL to obtain C and H using T . 

pw w w 

47. CALL SEVAL(l, TW, CPW , HW ) 

The enthalpy gradient along the wall: 

48. IIWP = (EJG)(CPW)(TWP) 

49. If TW ^ TAW, go to 54 

If TW > TAW, go to 50 

Write T and T . 
w aw 

50. WRITE TW, TAW 
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Write the axial distance x, Mach number M , momentum 

e 

thickness 0 , and energy thickness </> at the point where 

T exceeds T 
w aw 

51. WRITE X, ZME, THETA, PHI 

52. WRITE error indication message 
Stop the calculation by calling QUITS. 

53. CALL QUITS 

Save the wall enthalpy H . 

w 

54. A = HW 

Stagnation enthalpy minus wall enthalpy: 

55. B = HO - HW 

Save the minus sign of dynamic enthalpy ( ~U 2 /2) . 

G 

56. C = - UE202 

Save the free-stream temperature T . 

e 

57. TFINT = TE 

Call subroutine ZETAIT -to calculate the shape parameter 

1/h 

U = (A/6 ) ] and boundary layer thickness A, 6 , and 6* 

at point x, for given values of 0 and <p . 

58. CALL ZETAIT 

Save the value of p U /u . 

9 n 1 


59. CREY = (RHOUE) /(ZMU) 


Obtain the Reynolds number based on momentum 
thickness 0 . 


60. RTHE = (CREY) (THETA) 

The Reynolds number based on energy thickness <p : 

61. RPHI • (GREY) (PHI) 


Adiabatic 7' smperature T divided by free-s.ream 

aw 

temper.' : 

62. = (TAW)/(TE) 


Save (l - m) power of the above value 

63. ERASE 2 = (ERASEl)^ 1 * 0 " ZMVIS ) 


d-m 
T \ 
aw | 

T / ‘ 

e/ 


Set 

64. ERASE3 = ( 17. 2) (TO - TAW)/(TAW) 

65. ERASE4 = (305.0)(TE - T0)/(TAW) 

66. ICFCK = 0 


Save the guess value of skin friction coefficient 


fagt 


= ( 0 . 002 ). 


67. CFA = CFAGT 


Save the Reynolds number based on momentum thickness. 


68. RSUB = RTHE 

The following calculation, down to step 95, determines the 
friction coefficient by iterations. 
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69. Do 95, I ~ 1, 50 


Save 


70. CFAG - CFA 

Check the value of CFAGTP. 

71. If CFAGTP (ICFCII + l) = 0, goto 110 
If CFAGTP (ICFCII + l) * 0, go to 72 

Calculate the following expression from the relationship 
PM 

fr r \ = - e e (c \ tj 

v f e’ guess p JU ' fa' guess V 
aw aw 

72. CFRT = (ERASE2)(CFAG)(RSUB) 

The turbulent skin friction coefficient C* is obtained from 

subroutine CFEVAL, in which empirical relations between 
C^ and (C^Rg.) are tabulated. 

73. CFBAR = CFEVAL(CFRT) 

Sublayer temperature T divided by adiabatic wall 

temperature T : 
aw 

74. TSOTAW = 1.0 + (ERASE3) (SQRT({ CFBAR)/ (2.))) 

+ (ERASE4) ( (CFBAR) /(2.)) 

Check the sign of T /T 

s aw 

75. If TSOTAW > 0, go to 83 
If TSOTAW 0, go to 76 

Write error message. 

76. WRITE error message 
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Write X, M , 0, and <p in the case of BARPRO FAILURE. 

c 

77. WRITE X, ZME, THETA, PHI 

78. WRITE cause of error message 
Set 

79. CFAGTP ® (ICFCH + l) = 0.0 

80. CFAG =0.0 

81. CFA =0.0 

82. Go to 101 

Save the adiabatic skin friction coefficient 

83. CFA = (CFBAR)/(ERASEl)/(TSOTAW) ZMVIS 

Check the tolerance in the present C^ - C^R - 
iteration loop. 1 

84. If !(CFA-CFAG)/(CFAG)I < TOLCFA, go to 100 
If | (CFA-CFAG) /(CFAG) I S TOLCFA, go to 85 

Check the number of iteration. 

85. If I 2 2, go to 89 
If I < 2-, go to 86 

Set 

86. Z4 = CFA 

87. Z2 = CFAG 

88. Go to 95 

Save the following values. 
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89. 

Z3 = 

Z4 

90. 

Zl = 

Z2 

91. 

Z4 = 

CFA 

92. 

Z2 = 

CFAG 

93. 

Z'5 = 

(Z4 - Z3)/(Z2 - Zl) 

94. 

CFA 

= (Z4 - (Z5)(Z2))/(1. 


Go to 69 to continue the iteration. 

95. Continue 

Write that the skin friction coefficient could 
not be obtained. 

96. WRITE error message 

97. WRITE X, ZME, THETA, PHI 

98. WRITE Zl, Z2, CFA, Z3, Z4 

99. WRITE cause of error message 

Save C, . 
fa 

100. CFAG = CFA 

101. If ICFCH > 0, go to 111 
If ICFCH ^ 0, go to 102 

Save 

102. CFAGT = CFA 

103. CF = CFAG 

104. CFA = CFAGP 
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105. RSUB = RPHI 


SetICFCH equal one. 

106. ICFCH = 1 

Check whether the wall temperature is adiabatic or not. 

107. If ITWTAB £ o, go to 69 
If ITWTAB < 0, go to 108 

Set 

108. CH = 0.0 

109. Go to 113 

Check the value of ICFCH. 

110. If ICFCH ^ 0, go to 106 
If ICFCH > 0, go to 113 

Set 

111. CFAGP = CFAG 
Calculate the Stanton number. 


112. CH = ((PHl)/(THETA)) ZNSTAN ((CFAG)/(2.))/(l. - (5.) 
(SQRT( (CFAG)/(2. 0) ) (CHPARl)) 


Save 


1 d <» e u e> 

p U dx 
e e 


113. ERASE1 - (RHOUEP)/(RHOUE) 


Put 


+ 6 */e 


u 


dU 

e 

dx 


as 


49 


114. ERASE 2 = (UEP)(l.O + DKLSOT)/(UE) 

Call subroutine XNTERP to obtain the radius r and 
its gradient dr/dx at the axial distance x. 

115. CALL XNTERP (X, YR, YRP, IYX, XITAB, YITAB, IXTAB, 

CYX, IXPOS) 


Calculate the value 



116. DARC = SQRT(l.O + (YRP) 2 ) 
Put 


117. CDFORC = ( (RHOUE) /(G) ) (UE) /( DARC) (CF) /( 2. ) 

Check the geometry indicator EPSZ ( = 0: Two-dimensional 
planar flow, = 1: Axisymmetric flow ) . 


118. If EPSZ ^ 0.0, go to 120 
If EPSZ > 0.0, go to 119 


Save the value 


P U 
e e 


d (p e u e ) 

dx 


1 dr 
r dx 


119. ERASE 1 = ERASE 1 + ( EPSZ ) /(YR) ( YRP) 


The gradient of momentum thickness 

120. THE TAP = (CF)/;2.0)(DARC)- 

(THETA)(ERASE2 + ERASEl) 



Set 


121. ERASE2 = HO - HW 


The gradient of energy thickness 
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132. QDA = (ERASE1)(QW) 

133. DFORCE = (ERASEl) (CDi'ORC) 

134. DFLAT = 0.0 

135. Go to 139 
Set 


136. QDA = (QW)/(2.) 

137. DFORCE = (CDFORC)/(2. ) 

138. DFLAT = (DFORCE) ( Y RP) 

139. YOARC = Y2ARC 

140. Y2ARC = DARC 
Check the value of DX. 

141. If DX s o.O, return 


Call subroutine XNTERP to obtain the radius r = (ERASEll 
gradient | - (ERASE2) corresponding (q < ~ 

H2. CALL ~(X ^ £ 

XITAB. YITAB, KTABTcYxrSEOS) 
143. Y1ARC - SQRTfl.o + (ERASE2) 2 ) 

Increment of contour length: 

DXLARC = (DX)(yoaRC + (4.)(Y1AR C ) + Y2ARC)/(e. 
Contour length: 


145. XLARC = XLARC + DXLARC 
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rrr i -r — 


Integrated rate of heat transfer to the wall : 

146. SUMQDA = SUMQDA + (DXLARC) (QDA + QDAO) 

Drag in axial direction: 

147. FORCE = FORCE + (DXLARC) (DFORCE + DFORCO) 
Drag normal to the axial direction: 

148. FLAT = FLAT + (DXLARC) (DFLAT + DFLATO) 
Return to the main routine. 

149. Return 

Store. the results obtained in subroutine BARPRO for printout. 
Save the Reynolds number based on the contour length. 

150. RXLN = (CREY) (XLARC) 

Reynolds number based on the displacement thickness: 

151. RDLS = (CREY)(DELSTR) 

Check the value of the shape factor £ . 

152. If ZETA ^ 1.0, go to 160 
If ZETA < 1.0, go to 153 

Set 

153. 1=1 

Save the following integrals. 

154. Z1 = ZI4 


155. Z2 = ZI5 


156. Z3 = ZI6 


157. Z4 = ZI7 

158. Z5 = ZI1P 

159. Go to 166 
Set 

160. I = 6 

Save the following inegrals. 

161. Z1 = ZI1 

162. Z2 = ZI2 

163. Z3 = ZI3 

164. Z4 = ZI2P 

165. Z5 « ZI3P 
Print out 

166. WRITE heading for output values 

167. WRITE X, ZME, DELTA, HG, ZETA, CF 

168. WRITE XLARC, TE, BDELTA, QW, ZINTPR (i), Zl, CH 

169. WRITE YR, TW, DELSTR, SUMQDA, ZINTPR ( I + l), 

Z2, RTHE 

170. WRITE YRP, TAW, THETA, FORCE, ZINTRP (I + 2), 

Z3, RXLN 

171. WRITE ZMEP, PHI, FLAT, ZINTPR (I + 3), Z4, RPHI 

172. WRITE UE, DELSOT, ZINTPR (I + 4), Z5, RDLS 
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173. WHITE PE 


174. Return 

Table 4 gives subroutine BAltPliO nomenclature. 


TABLE 4. SUBROUTINE BA IIP 110 NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

A 

Saved value of wall 
enthalpy II 

ft 2 /sec 2 

/SAVED/, 54 

B 

Stagnation enthalpy minus 

wall enthalpy H - II 
o w 

ft 2 /sec 2 

/SAVED/, 55 

BDELTA 

Temperature thickness A 

ft 

/OUTPUT/, 

168 

C 

Minus sign of dynamic 
U e 

enthalpy- — 

ft 2 /sec 2 

/SAVED/, 56 

CDFORC 

Local drag force per 
unit area 

lbf/ft 2 

117, 133, 137 

CF 

Skin friction coefficient 

' 111 " 

/OUTPUT/, 
103, 117, 120, 
167 

CFA 

Adiabatic skin friction 
coefficient C 

fa 


67, 70, 81, 83, 
84, 86, 94, 

100, 102, 104 

CFAG 

Guess value of adiabatic 

skin friction coefficier C„ 

fag 


70, 72, 80, 84, 
87, 100, 103, 
111, 112 

CFAGP 

Skin friction coefficient 
obtained by using the 

Reynolds number R , 

0 


/INTER/, 104, 
111 

CFAGT 

Initial value of skin 
friction coefficient 


/INTER/, 
EQUIV, 67, 
102 

CFAGTP 

Array equivalenced to 
CFAGT and CFAGP 


DIM, EQUIV, 
71, 79 

CFBAR 

Skin friction coefficient 
obtained from empirical 
relationship 


73, 74, 83 



TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (C -ntinued) 


Symbol 


Description 


Units 


Reference 


CERT 


CHPAR1 


CPW 


CREY 


CTWX 


DARC 


DELSTR 


DELSOT 


Turbulent skin friction 
coefficient multiplied 
by experimental Reynolds 
Number Cj,R- 


Stanton number 


- P r + In 

\ 5Pr + 1 


Array of parabola 
coefficients for Mach 
number table 


Specific heat in free 
stream C 

pe 


Specific heat at the 
wall temperature C 

pw 


Mass flow density 

divided by viscosity 

in free- stream 

pU/ft 
e e e 


Array of parabola 
coefficients for the 
wall temperature table 


Array of parabola 
coefficients for the 
free stream temperature 
table 


Array of parabola 
coefficients for the nozzle 
r a dius table 


1 + (dr/dx) 2 


Displacement thickness 
6 * 


Displacement thickness 
divided by momentum 
thickness 6 */0 


72, 7.1 


/OUTPUT/, 
108, 112, 122, 
125, 168 


/INTER/, 112 


/LOOKUP/ 2, 
20, 28 


Btu/( lbm - °R) 6, 7, 8 


Btu/( lbm 


47, 48 


59, 60, 61, 
150, 151 


/LOOKUP/, 

46 


LOOKUP/, 4 


/LOOKUP/, 
115, 142 


116, 117, 120, 
122, 140 


/OUTPUT/, 
151, 169 


/OUTPUT/, 
114, 172 










































TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

DELTA 

Velocity thickness 6 

ft 

/OUTPUT/, 167 

DFLAT 

One-half the local di’ag 
force normal to x-axis 

ibf/ft 2 

128, 130, 138, 
148 

DFLATO 

Saved value of DFLAT 

lbf/ft 2 

129, 148 

DFORCE 

One-half the local drag 
force (for axisymmetric 
flow) or 

(for two-dimensional 
planar flow) 

lb!7ft 

lbf/ft 2 

133, 147 

128, 137, 138, 147 

DFORCO 

Saved value of DFORCE 

lbf/ft 2 

128, 147 

DX 

Weighted difference of 
table values of axial 
distance 

ft 

/INTER/, 141, 
142, 144 

DXLARC 

Increment of the length 
along contour 

ft 

144, 145, 146, 147, 
148 

DXRHO 

One-tenth the difference 
between axial distance 

ft 

/INTER/, 13_, 28, 
31 

EPSZ 

Geometry indicator 
EPSZ = 0: two- 

dimensional 
planar flow 

EPSZ = Is axisymmetric 
flow 


/INPUT/, 118, 
119, 130 

ERASE 1 

1. Free -stream Mach 
number M 

e 

2. Saved value of T / T 

3. Saved -value of aw e 

1 d <»e U e> 

p U dx ° r 

e e 

i d (P U ) . . 

1 eel dr 

p U dx r dx 

e e 

4. Trr 

5. Contour radius r 

ft 

ft 

20, 21, 28, 29 

62, 63, 72, 83 
113, 119, 120. 122 

131, 132, 133 
142 


57 



TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 


ERASE 3 


Description 


ERASE2 1. Gradient of free - 

stream Mach number 

dM /dx . 
e 

2. Saved value of 


3. Saved value of 


1+6 */Q dU 


4. Enthalpy difference 

H - H 
o w 

5. dr/dx 


1. Specific heat at 
adiabatic wall 
temperature C 

paw 

2. Saved value of 17.2 

(T - T )/T 
o aw aw 



Reference 


20, 28 


114, 120 


121, 122 
142, 143 


36 


64, 74 


ERASE4 



Saved value of 305 

(T - T )/T 

e o aw 

Conversion factor 

between thermal and 

work units 

Conversion factor 

between thermal and 

work units multiplied by 

acceleration of gravity 

used as a proportionality 

constant gJ 


(ft-lbf)/Btu 


65, 74 


''INPUT/, 125 


lbm/( BUi-sec 2 ) /C SE VA L/, 7 , 
48 
















TABLE 4. SUBROUTINE BARPRO NOMENCLATURE ( Cont inued) 


Symbol 

Description 

Units 

Reference 

FLAT 

Force normal to x- 
direction_for two- 
dimensional planar flow 

Ibf 

/output/, 

148, 171 

FORCE 

Drag force in axial or 
x-direetion 

ibf 

/output 77 

147, 170 

G 

Acceleration of gravity 
used as a proportionality 
constant 

lbm7lbf 
ft/ sec 2 

/input/, 117, 

125 

GAME 

Specific heat ratio in 
free stream 

— 

8 

HO 

Stagnation enthalpy II 

0 

ft 2 /sec 2 

/CSEVAL/, 9, 
55, 121 

HAW 

Adiabatic wall 
enthalpy H 

aw 

ft 2 /sec 2 

.‘35, 56, 39, 
122, 125 

HE 

Enthalpy of free stream 
in work units H 

e 

ft 2 /sec 2 

/INTER/, 6, 
9, 35 

HEP 

Enthalpy gradient 

dH /dx 
e 

ft/sec 2 

7, 11, 40 

HG 

Heat transfer 
coefficient h 

g 

Btu/(ft 2 -sec° R) 

/OUTPUT/, 
126, 167 

HW 

Enthalpy at the wall H 

w 

ft 2 /sec 2 

/INTER/, 39, 
43, 47, 54, 55, 
121, 122, 125 

HWP 

Enthalpy gradient 

dH /dx 
w 

ft/sec 2 

40, 44, 48 

I 

Subscript counter 

— 

153, 160, 168- 
172 

ICFCH 

Option indicator 

— 

66, 7.1, 79, 
101, 106, 110 

IMX 

Mach number table entry 
indicator and saved 
subscript counter 


/LOOKUP/, 2, 
20, 28 

IND 

Program loop control 
indicator 

~-- J — 

CALL, 1, 123 


59 



















TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 




Symbol 

Description 

U nits 

Reference 

ITV7TAB 

= - Is T = adiabatic 
w 

wall temperature 

= 0: T - constant . 

w 

= 1 5 IXTAB values of 

T will be input 
w 

— 

/INPUT/ 37, 
107, 124 

ITWX 

Wall temperature table 
entry indicator, and-saved 
subscript counter . 


/LOOKUP/, 46 

ITX 

Free -stream temperature 
table entry indicator and 
saved subscript counter 


/LOOKUP/, 4 

IXPOS 

Array position indicator 

— 

/LOOKUP/, 2, 
4, 20, 28, 46, 
1 . 15 , 142 


IYX 


PE 


PHI 


PPHIP 


PIE 


PRE103 


QDA 


QDAO 


60 


x, y, and Mach number 
tables (4 ^ IXTAB 
=s 500) 


Nozzle radius table entry 
indicator and saved 
subscript counter 


Static pressure of the 
free stream 


Energy thickness.^ 


Gradient of energy 
thickness d£ /dx 


Circumferential 
constant 7r 


Cubic root of Prandtl 
number 


One -ha If the local heat 
transfer to the wall 

( axisymmetric flow) 
( two-dimensional 
planar flow) 


Saved value of QDA 


lbf/ft* 


ft 


Btu/(ft-sec) 

Btu/(ft 2 _sec) 


Btu/(ft-sec) 


20, 28, 46, 
115, 142 


/LOOKUP/, 
115, 142 


/OUTPUT/, 5, 
12, 173 


/OUTPUT/, 

51, 61, 77, 97, 
112, 122, 171 
/INTER/, 122 


/INPUT/, 131 


/INTER/, 35, 
40 


127, 132, 146 
127, 136, 146 


127, 146 









TABLK 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

QW 

Local heat transfer rate 
to wall 

Btu/(ft 2 -sec) 

/OUTPUT/, 
125, 126, 132, 
136, 168 

RBAR 

Gas constant in work 
units 


/INPUT/, 12, 
30 

RJ)LS 

Reynolds number based 

on displacement thickness 

R„ - P U 6*/M 
e e e 


151, 172 

RIIOE 

Free -stream density p^ 

lb m/ft 3 

/INTER/, 12, 
32, 33 

RIIOEP 

Density gradient dp /dx 

lb m/ft 4 

14, 31, 33 

RIIOUE 

Mass flow density p^U^ 

lbm/(ft 2 -sec) 

/INTER/, 32, 
59, 113, 117 

RIIOUEP 

Gradient of mass flow 
density d(p^U^) /dx 

lbm/(ft ,3 -sec) 

33, 113 

ROJ 

Gas constant in thermal 
units 

Btu/ ( lbm-° R) 

/CSEVAL/, 8 

RPIII 

Reynolds number based 

on energy thickness 

R . = P U ■p/p. 
p c e e 


61, 105, 171 

RSUB 

Saved value of Reynolds 
number R^ 

Saved value of Reynolds 

number R , 

0 


68, 72, 105 

RTHE 

Reynolds number based 

on momentum thickness 

R n = P U 0/p 
0 c c e 


60, 68, 169 

RXLN 

Reynolds number based 

on arc length 

R t = p U h/p 
L e e e 


150, 170 

SUMQDA 

Integrated heat transfer 
rate to the wall up to 
station x 

(axisymmetric flow) 
( two-dimensional 
planar flow) 

Btu/sec 
Btu/( sec -ft) 

/OUTPUT/, 
149, 169 


61 








TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 


Description 


Units 


Reference 


TO 

St; i gnu t i on to mpc ra tur e 
T 

o 

° R 

/INPUT/, 34 

TAW 

Adiabatic wall tempera- 
ture T 

aw 

° R 

36, 38, 49, 50, 
62, 126, 170 

TE 

Frcc-stx’eam tempera- 
ture T 

e 

° R 

/OUTPUT/, 4, 
5, 6, 12, 34, 
57, 62, 168 

TEP 

Temperature gradient in 

free -stream dT /dx 
e 

0 R/ft 

4, 7 

TFINT 

Saved value of free- 
stream temperature 

°R 

/COFIIF/, 57 

THETA 

Momentum thickness 0 

ft 

“/output/, 


THETAP 

Gradient of momentum 
thickness d0/dx 

TITAB 

Free ..stream tempera- 
ture table related to 
IXTAB array. This 
table is determined by 
subroutine GETPT via 
BARSET. 

tolcfa 

Tolerance in C" f - C f R- 
iteration 

TSOTAW 

Sublayer temperature 

divided by adiabatic wall 

temperature T /T 

s aw 

TW 

Wall temperature 

TWP 

Temperature gradient on 
the wall 

TWTAB 

Wall temperature table 
related to IXTAB array 


51, 60, 77, 

97, 112, 120, 

170 

/INTER/, 120 

/TABLES/, 4 



/INPUT/, 84 


74, 75, 83 


T 

42, 46 
49. 50 


46, 48 


TABLES/, 42, 
43, 46 














TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

UE 

Free -stream velocity 

ft/sec 

/OUTPUT/, 
10, 11, 32, 
33, 10, 114, 
117, 172 

UE202 

Square of free -stream 

velocity divided by two 

(dynamic enthalpy) 

U 2 /2 
e 

ftVsec 2 

9, 10, 56 

UEP 

Velocity gradient along 
free stream 

l/sec 

11, 33, 40, 114 

X 

Axial distance or 
distance in free-stream 
direction 

ft 

/OUTPUT/, 

2, 4, 16, 20, 
24, 28, 46, 
51, 77, 97, 
115, 142, 167 

XIBASE 

First value in axial 
distance table 

ft 

./INTER/, 16 

XIEND 

Last value in axial 
distance table 

ft 

/INTER/, 24 

XITAB 

Table of EXTAB values 
(axial distance x) in 
monotonically increasing 
order 

ft 

/TABLES/, 2, 
20, 28, 46, 
115, 142 

XLARC 

Arc length of contour 
corresponding to x 

ft 

/OUTPUT/, 
145, 150, 168 

YOARC 




139, 144 

Y1ARC 

7* + (U ' 

corresponding to 
(-?) 


143, 144 

Y2ARC 


— 

139, 140, 144 

YITAB 

Nozzle radius or contour 
height table related to 
IXTAB array 

ft 

/tables/, " 

115, 142 

YR 

Radius or height of 
contour r 

ft 

/OUTPUT/, 
115, 119, 
131, 169 



















































TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

YRP 

dr/dx 

— 

115, 116, 119, 
138, 170 

Z1 

1. Saved value of free- 
st ream density p 

e 

mm 

17, 22, 31 


2. Saved value of C. 

fag 

3. Saved value of 
I t = (ZIl) or 
I 4 = (ZI4) 


90, 98 


■ 

/OUTPUT/, 
154, 161, 168 

Z2 

1. Saved value of p 

e 

lb m/ft 3 

25, 30, 31 


2. Saved value of 
C tag (CFAG) 

— 

87, 92, 94, 98 


3. Saved value of integral 
I, = (ZI2) or 
I 5 = (ZI5) 


/OUTPUT/, 
155, 162, 169 

Z3 

1. Saved value of one or 
point five 

— — 

18, 23, 26, 31 


2. Saved value of C^ a 

— 

89, 98 


3. Saved value of integral 
I 3 = (ZI3) or 
I 6 = (ZI6) 

■■■■ 

/OUTPUT/, 
156, 163, 170 

Z4 

1- Free-stream pressure 
P 

e 

lbf/ft* 

21, 22, 29, 30 


2. Saved value of 
C fa (CFA) 

— 

86, 91, 94, 98 


3. Saved value of integral 
U = (ZI2P) 

— 

/OUTPUT/, 
157, 164, 172 

Z5 

. 1. Free-stream tem- 
perature T 

e 

°R 

21, 22, 29, 30 


2. Saved value of 

(Z4 - Z3)/(Z2 - Zl) 

— 

93, 94 


3. Saved value of integral 
I\ = (ZI1P) or 
IJ = (ZI3P) 


/OUTPUT/, 
158, 165, 172 

ZETA 

Shape factor: 
t = (A/6) n 


70UTPUT 77 

152, 167 



























TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

ZI1 

Integral 

h - / J-S n (l - S)dS 
o 


/OUTPUT/, 161 

ZI2 

Integral 

h - S f- s"dS 

o 

1 

/OUTPUT/, 162 

ZI3 

Integral 

1 3 = / — S n_1 dS 
o p e 

' 

/OUTPUT/, 163 

ZI4 

Integral 

>4- / T-S n (l - S)dS 
o 


/OUTPUT/, 154 

ZI5 

Integral 

i 5 = / ~s n (i - S)dS 


/OUTPUT/, 155 

ZI 6 

Integral 

Ie = / — S n dS 
o p e 


/OUTPUT/, 156 

ZI7 

Integral 

I, - / s n ds 

£ P e 


/OUTPUT/, 157 

< 

ZI1P 

Integral 

1 

i’i = / 7 - w n (i - W)dW 
0 p e 


/OUTPUT/, 158 































TABLE 4. SUBROUTINE BARPRO NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

ZI2P 

Integral 

1/t 

I } = f — W U (1 - W)dW 
o P e 


/OUTPUT/, 

164 

ZI3P 

Integral 

1 

x 3 = / ~ W n-1 (l - W)dW 

l/£ P e 


/OUTPUT/, 

165 

ZINTPR 

Array of Hollerith headings 
for printout of integrals 


DATA, 
DIM, 168, 
169,’ 170, 
171, 172 

ZME 

Mach number of free-stream 
M 

e 


/OUTPUT/, 
2, 4, 51, 

77, 97, 167 

ZMEP 

dM /dx 
e 

1/ft ... 

2, 171 

ZMTAB 

Mach number table related 
to EXTAB array 

— 

I 

ZMU 

Viscosity of free stream 
p e 


34, 59 

ZMUO 

Stagnation viscosity n 

0 

lbm/(sec-ft) 

/INPUT/, 

34 

ZMVIS 

Exponent in viscosity - 
temperature law 



ZNSTAN 

Interaction exponent in 
Stanton number relation S' 


/INPUT/, 

112 
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SUBROUTINE BARSET 


Subroutine BARSET computes program constants and sets up specific 
heat C and enthalpy H versus temperature T tables and tables of pressure 

P and temperature T as a function of axial distance x for the inviscid flow. 
The entropy S for T and P can be determined. 

Evaluation of enthalpy is as follows: 


II(T) 


H(T) 


H(T) 


X 

“ Hj + / C (t)dt , 


T. 


T. 


* A 

= / C (t)dt + J' C (t)dt , 


o P T. P 

1 


T 1 T 2 % 

= / C (t)dt + f C (t)dt + y C (t)dt + 

^ P rp P rp P 


T 

+ f c (t)dt . 

T. p 
1 



The first integral is evaluated assuming that the 
linear from T = 0 to T = T t . 


C 

P 


- T curve is 


The value of C at T = T< is C . . The value of dC /dT at T = T. 

P 1 pl P 1 

is BCPj. Thus for 


0 < T < Tj , 


C 

P 



o 


= C + BCP^T - T x ) , 
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Therefore, one obtains 


1 BCP 

/ C (t)dt = C T — i t 2 

J Q P Pi 1 2 1 


The remaining integrals can be evaluated gradually from the expression for C 



+ DCP.(t - T.) 3 ]dt , 


where 



Therefore, 
T 


i+1 


(Tu, - T -) 2 


/ C (t)dt = C (T - T ) + ™ ' i+1 *' 


T. 


pi'* i+1 ~ *i' 1 BCP i 2 


(T i+l - T i ) 3 (T i+1 - T ) 4 

+ CCP. — — i- + DCP — — i— 

13 i 4 * 


and 


T 

/ c p (t)dt - C pi (T -T,) + BCP. 2 
(T - T^ 4 


(T - T.) 2 


r (T - T ) 

+ CCP 1 


3 


i 3 


A • 


+ DCP, 



By assuming that 



H. 

1 




C (t)dt 
P 


» 


then. 

(T - T.) 2 

H(T) = H. + C pi (T - T.) + BCP. 5-^- 


(T-T.) 3 (T - 1 

+ ccp i— H- + DCP i— - 


The values of H are calculated and stored in BARSET under the name 
i 

HTAB (I). 


Evaluation of entropy is as follows. 



- Rln (t) • 

Arbitrarily S = 0 can be set at T = T t and P = P Q where Tj is the' first 
point of the C - T table and P is the stagnation pressure. 



dt 


Tj C (t) 

-8 + - 4 - 

o £ t 
o 


Thus the entropy will be evaluated from the above equation. 


T C (t) 


s ’l~ t 


dt - R In 


Ti 

By definition 



then 


P’ R’ = R In ( 1 

iP o> 


T C (t) 

S = / JL- dt _ P <R. . 

x, t 


T 2 C (t) T, C (t) 

S - / J^dt+ / j£-dt + ... 

Ti T 2 


T 

I 

T 


i C (t) T C (t) 

+ J Jt— dt + J J£— dt - P’R’ 

m ^ m t 


i -1 


i 


T 


i-1 j+1 C (t) 

S = Yj J - 7 — dt + f -X — dt - P» R* 

l 


j = l T. 


T C (t). 
" £_ 
t 


1 


i 


From the expression for C^, 


i-1 T i + 1 C + BCP (t -T.) + CCP.(t - T.) 2 + l)CP.(t -T.) 3 

s = y; f ' -El 3 1 —i j j — i- dt 

j=l T 1 

J j 


T C + BCP (t - T ) + CCP (t - T.) 2 + DCP.(t - T.) s 

/» Y\ 4 i ' S' S ' 4 7 1 ' 1 7 


I 

T. 


Pi V r i i i 


— dt - P' R’ . 


T 

i-1 1 + 1 

J 

j=l T. 


C BCP (t - T ) CCP (t - T.) 2 DCP.(t - T.) 3 

_PJ + 3 L + 3 L_+ 2 2 — 

t t t t 


dt 


J 


T 

♦ J 

T. 

l 


C . BCP.(t - T.) CC?.(t - T.) 2 DCP.(t - T.) 3 

_£} + i 1 + 1 1 + - i- 

t t t t 


dt - P’R* 


i-1 

s=E 

j-i 


c lnt + BCP (t - T. In t) + CCP. 
PJ J J j 


(I - 2tT , 


+ T 2 In t j 


+ DCP 


(?- 


3t 2 T. ^ 

— 2~ 3 + 3cT 2 - T 3 In t j 


T. 


j+1 


+ jc p . lnt + BCP.(t - T.lnt) + CCP.^| - 2tT. + T 2 


4.3 3t 2 T 

+ DCP.I — + — rr- + 3tT 2 - T 3 In t 
i\ 3 2 l i 


- P* R* 
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i-1 


: .In T + BCP,(T , - T In T ' 
pi j+i y j+i j j-i-r 


r p2 

+ CCP l-^ - 2T, T -I- T 2 In T 

;i\ 2 .1+1 .1 j J+Jy 

rT 2 3T 2 T. 

+ DCP.l — + 3T../J' 2 - T* In T. 


1\ 3 


} 


.1+1 i j jii 


i’?. 


- c p j T j - BCP j (T j “ T j T i ) - CCP ilT - 2 T j T j + T - ln 

/ T 3 3T 2 T \‘ 

- DCP j\-jf — jr 1 + st.t 2 - t 2 hi t J 

+ C . In T + BCP.(T - T. ln T) + CCP Z-- 2 - - 2T-T + T 2 ln 
P 1 l i . i\ 2 i i ) 

/V 

i\3 


3T 2 T 


+ DCP 


-r — + 3TT 2 - T? ln T , 

2 l l ) 


/t 2 

- C p .ln T. - BCP i (T i - T. In T.) - CCP.^-~ - 2T.T^ + T 2 ln T 

^>3 rj>2rp 


i\ 3 


1 l 

2 


- DCP.l -rr - 3 + 3T.T 2 - T 2 In T. 


l l 


- P’R’ 




then 


S : V 


i-1 
S' 
_/ 
i l 


In j' 1 


BAHBi. + (T. 4 - T.) BARB2. 
1. J J+l J J 

J 


DCP. 

i (T 2 4 - T 2 ) BARB3. -i- (T 2 . T 3 ) — — l 
J+l J J j+l J 


DCP. 

+ In T BARB1, + T BARB2 . + T 2 BARB3. + T 3 — — i 
l l l 3 


DCP. 

- In T. BAHBI. - T. BA11B2. - T? BARB3 . - T? — — ^ 
l li li l l 3 


The terms G1 and G2 are defined as follows: 


i-i 

j=i 


In BARB1 . + (T. - T.) BARB2. 

T. J J+l J J 


DCP. 

+ (T 2 - T 2 ) BARB3 . + (T 3 - T 3 ) — — i 

J+l J J J+l J 3 


DCP. 

G2 = In T. BARB1. + T. BARB2. + T? BARB3. + T? — — i 
l li li l i 3 


Then 

T 3 DCP. 

S = Gi + In T BA RBI . + T BARB2. + T 2 BARB3. + - — - 

i l l 3 


- G2 


By letting 

G. = Gl - G2 , 

l ’ 


P'R' . 


P'R' . 
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ono finally obtains 


S = G. + BARB1. In T + BARB2.T + BARB3.T 2 + DCP, ?! - P’ R' . 
1 l 1 i i 3 

Tho terms Gl, G2, G,, BARB1., BARB2., BARB3. are evaluated in 

l i i i 

BARSET. 

COMMON BLOCKS 

COMMON blocks CSEVAL, INPUT, and TABLES are used. 

TBL SUBROUTINES 

Subroutine DIRECT calls BARSET. 

BARSET calls subroutines BMFITS, SEVAL, QUITS, and GETPT. 
FORTRAN SYSTEM ROUTINES 

FORTRAN library routine ALOG is used. 

CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL BARSET 
SOLUTION METHOD 

Set the circular constant it. 

1. PIE = 3.14159265 

Multiply sea-level acceleration with the conversion factor 
relating thermal and work units. 

2. FJG = (FJ)(G) 

Divide the specific gas constant by the conversion factor 
relating thermal and work units. 

3. ROJ = (RBAR)/(FJ) 


Set TMAX equal to the free-stream stagnation temperature 


4. TMAX = TO 


Set indicator II for nominal entry. 

5. II = l 

Test whether adiabatic wall temperature (ITWTAB = -l), 
constant wall temperature (ITWTAB = 0 ) , or tabular wall 
temperature (ITWTAB = l) is to be used. 

6. If ITWTAB > 0, go to 7 
If ITWTAB = 0, go to 8 
If ITWTAB < 0, go to 12 

Set II equal to the number of points in x, y, and M tables (IXTAB). 

7. II = IXTAB 

Check and save the maximum value of the wall temperature tabulated 
(TWTAB). 

8. Do 11 I = 1, II 

9. If TWTAB (I) < TMAX, go to 11 
If TWTAB (I) * TMAX, go to 10 

Save the maximum value of TWTAB. 

10. TMAX = TWTAB (i) 

11. Continue. 

Test whether constant specific heat calculation (ICTAB = 0 ) 
or specific heat table is required (ICTAB > 0 ) . 

12. If ICTAB = 0, go to 47 
If ICTAB * 0, go to 13 


Check whether the temperature input value (TMAX) exceeds the table 
upper limit. 

13. If TMAX < TCTAB (ICTAB), go to 16 

If TMAX 2= TCTAB (ICTAB) , go to 14 

WRITE TCTAB (ICTAB) and TMAX, when the table upper limit 
TCTAB (ICTAB) is exceeded. 

14. WRITE TMAX, TCTAB (ICTAB) 

Stop the calculation by calling QUITS, 

15. Go to 49 

Save the number of points in C^ versus T table (ICTAB). 

16. NOCTAB = ICTAB 

The following calculation, down, to step 43, is only used for specific 
heat polynomial, enthalpy, and entropy equations. Call subroutine 
BMTAB to determine polynomial coefficients BCP, CCP, and DCP. 

17. CALL BMFITS (TCTAB, CPTAB, ICTAB, BCP, CCP , DCP ) 
Calculate temperature and various powers thereof. 

18. TIE1 = TCTAB (l) 

19. TIE 2 = (TCTAB(l)) 2 
-20. TIE3 = (TCTAB(l)) 3 

Calculate the enthalpy at Tj. 

21. HTAB(l) = (CPTAB(l)) (TIEl) - (BCP (l) ) (TIE2)/(2. 0) 

A term in entropy equation 

22. BARBl(l) = CPTAB(l) - (BCP(l) ) (TIEl) + (CCP(l) ) (TIE2)- 

(DCP(1))(TIE3) 
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First derivative of the above equation with negative sign 


23. BARB2(1) = BCP(l) _ (2. 0) (CCP( 1 ) ) (TIEl) 
+ (3,o) (l)CP(l )) (TIIO2) 


Second derivative of the equation 22 times ( -0.25). 

24. BARB 3 ( 1 ) " ((CCP(l)) - (.3. 0) (DCP (l ) ) (TIEl)/(2. 0) 
First term of the entropy equation 

25. GTAB(l) = - (BARB 1 ( 1 ) ) ( Ln ( TIE 1 ) ) - (BARB2(l))(TIEl)~ 

(BARB3(l ))(TIE2) - (UCP(l ) ) (TIE3)/(3. 0) 

Set entropy summation term to zero. 

26. G1 = 0.0 

Store the enthalpy and coefficients of entropy equation at each 
temperature corresponding to - T table. 

27. Do 43, I = 2, ICTAB 

Save 71 = TIEl, Tl 2 = TIE2, Tl 3 = TIE3. 

28. TME1 = TIEl 

29. TME2 = TIE2 

30. TME3 = TIE3 

Set new values from C - T table. 

P 

31. TIEl = TCTAB (i) 

32. TIE2 = (TIEl) (TIEl) 

33. TIE 3 = (TIE1)(TIE2) 

Determine the temperature difference (T^ - T^ ^). 

34. DELT = TIE1-TME1 

Enthalpy at temperature T from C - T table. 

i p 
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35. IITAB(I) = HTAB(I - l) + (CPTAB(l - l))(DELT) 
+ (BCP (i _ l) ) (delt) 2 / ( 2 . 0 ) 

+ (CCP(l - 1)) (DELT) 3 /(3.0) 

+ (DCP (i - 1)) (DELT)V(4.0) 


Check whether the value of I equals or exceeds ICTAB. 

36. If I ^ ICTAB, go to 43 
If I < ICTAB, go to 37 


Determine the coefficients used in the entropy equation: 


S = GTAB(I) + BARBl(l) Ln T + BARB 2 (i) T + BARB3 (i) T 2 
+ DCP(l) T 3 /3 - R Ln(P/Po) . 

37. BARBl(l) = CPTAB(l) - (BCP (i) ) (TIEl) + (CCP(l)) (TIE2) 
- (DCP(l))(TIE3) 


First derivative of the previous equation 

38. BARB2(I) = BCP(l) - (CCP(l)) (TIEl)/(0.5) 
+ (3.0) (DCP (i) ) (TIE2) 


Second derivative of 37 times (-0.25)' 

39. BARB3(I) = (CCP(l) - (3.0) (DCP (i)) (TIEl)) /(2.0) 


Sum up terms in entropy equation. 

40. G1 = G1 + (BARB l(l - l))(Ln((TIEl)/(TMEl))) 

+ (BARB2(I - l))(DELT) 

+ (BARB3 (i - l) ) (TIE2 - TME2) 

+ (DCP (I - l) ) (TIE3 - TME3)/(3.0) 

41. G2 = (BARBl(l) ) (Ln(TIEl) ) + (BARB2(l) ) (TIEl) 

+ (BARB3(I))(TIE2) + (DCP(l))(TIE3)/(3.0) 

42. GTAB(I) = G1 - G2 


43. Continue 


Determine the number of input intervals minus one. 

44. IS = ICTAB - 1 
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are obtained from 


Stagnation enthalpy II ^ and specific heat 
subroutine SEVAL for a given stagnation temperature. 

45. CALL SEVAL( 1, TO, CPO , IIP) 

Calculation of specific heat ratio at stagnation point 



= C /(C 
po po 


- R/J) 


46. GAMO = (CPO) /(CPO - ROJ) 


When this specific heat ratio y^ is less than 1, check for inconsistant 
units and stop the calculation. Otherwise continue. 


47. If GAMO >1.0, go to 50 
If GAMO ^ 1.0, go to 48 


48. WRITE GAMO 


Subroutine QUITS stops the computation. 

49. CALL QUITS 

Set (y - l)/2 as 
o 

50. GM102 * (GAMO - 1. )/(2. ) 

51. GOGM1 = (GAMO) /(GAMO - 1.0) 

Save the stagnation pressure. 

52. POMAX = PO 

Generate the table of free-stream pressure P and temperature T^. 
P g and T^ are obtained for a given by using subroutine GETPT. 


53. Do 54, I = 1, IXTAB 


54. CALL GETPT(ZMTAB(I) , PITAB(l) , TITAB(T)) 


Calculation of the free-stream enthalpy in the case of C = constant. 
This calculation continues to step 85. p 

55. If ICTAB > 0, so to 85 

If ICTAB =£ 0, go to 56 


Set 


56. NOCTAB = 6 

Save NOCTAB minus one. 

57. NOCTM1 = NOCTAB - 1 
Save NOCTM1. 


58. IS = NOCTM1 


Compute the specific heat for stagnation condition C 

po 


59. CPO = (GOGMl)(RBAR)/(Fj) 


y R 

J(y 0 °" o 


Y g R 

Calculate the specific heat in work units C' = 

p y o - 1 

60. CJG ■ (CPO) (FJG) 


Compute stagnation enthalpy in wor k units. 

61. HO = (CJG) (TO) 

If the free-stream temperature exceeds the stagnation temperature, 
set the former as TMAX. 


62. Do 66, I = 1, IXTAB 

Set the free-stream temperature equal to the table value T . 


63. TE = TITAB(I) 

64. If TE s TMAX, go to 66 

If TE > TMAX, go to 65 

Save T . 
e 

65. TMAX = TE 

66. Continue 

Set the maximum temperature in the table to TMAX plus a hundred. 

67. TCTAB(NOCTAB) = TMAX + 100.0 
Set the first temperature table value to 10 -10 . 

68. TCTAB(l) = 1.0E-10 
Obtain NOCTM1 as a real number. 

69. Z1 = NOCTM1 

Definition of DELT = ™ AX + 100 - ( 1 »0 E - 1 0) 

5 

70. DELT = ( TCTAB(NOCTAB) - TCTAB(l) )/(Zl) 

Compute first term in entropy equation. 

71. ERASE 1 = -(CP0)(Ln(TCTAB(l)) 

Compute the coefficients for the specific heat polynomial 
and entropy equation. 

72. Do 84, I = 1, NOCTAB 
Save ERASE 1. 

73. GTAB(I) = ERASE1 


Set the specific heat C in C table to a constant value. 

po p 

74. CPTAB(I) = CPO 

Set coefficients in the C = polynomial equal to zero. 

P 

75. BCP(I) =0.0 

76. CCP(I) =0.0 

77. DCP(I) =0.0 

Set coefficients in the entropy polynomial. 

78. BARBl(l) = CPO 

79. BARB2(I) =0.0 

80. BARB3(I) =0.0 

Compute table values of enthalpy corresponding to TCTAB, 

81. HTAB(I) = (CP0)(TCTAB(I) - TCTAB (l) ) 

82. If I a NOCTM1, go to 84 
If I < NCOTM1, go to 83 

Generate the temperature table TCTAB. 

83. TCTAB (I + 1) = TCTAB (l) + DELT 

84. Continue 

Check whether variable wall temperature option (ITWTAB # 0) is used. 

85. If ITWTAB * 0, go to 87 
If ITWTAB = 0, go to 86 

Obtain the specific heat and enthalpy for constant wall temperature by 
calling the subroutine SEVAL. 

86. CALL SEVAL( 1, TWTAB(l), ERASE 1, TWTAB(2)) 
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Obtain the stagnation entropy from stagnation temperature and pressure 
by ealling the subroutine SEVAL. 

87. CALL SEVAL (0,T0, PO, SO) 


88. Return 

Table 5 gives subroutine BARSET nomenclature. 

TABLE 5. SUBROUTINE BARSET NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

BARB1 

Coefficient in polynominal 
equation 

Btu/( lbm-° R) 

/C SEVAL/, 
22, 25, 37, 
40, 41, 78 

BARB2 

First derivative of BARB1 
(negative sign) 

Btu/( lbm-° R 2 ) 

/CSEVAL/, 
23, 25, 38, 
40, 41, 79 

BARB3 

Second derivative of 
BARB1 times -0.25 

Btu/(lbm-°R 3 ) 

/CSEVAL/, 
24, 25, 39- 
41, 80 

BCP 

Coefficient in the C -T 
P 

relationship determined by 
BMFITS 

Btu/(lbm-°R 2 ) 

/CSEVAL/, 
17, 21-23, 
35, 37, 38, 
75 

CCP 

Coefficient in the C -T 
P 

relationship determined by 
BMFITS 

3tu/ ( lbm- 0 R 3 ) 

/CSEVAL/, 
17, 21-24, 
35, 37-39, 
76 

CJG 

Specific heat at stagnation 
condition in work units 

ft 2 /(sec 2 - 0 R) 

/CSEVAL/, 
60, 61 

CPO 

Specific heat at constant 
pressure at stagnation 
condition 

Btu/( lbm-° R) 

/CSEVAL/, 
45, 46, 59, 
60, 71, 74, 
81 

CPTAB 

Array of C values corre- 
P 

sponding to the values in 
the temperature table 

Btu/(lbm-° R) 

/CSEVAL/, 
17, 21, 22, 
35, 37, 74 

DCP 

Coefficient in the C -T 
P 

relationship determined by 
BMFITS 

Btu/( lbm-° R 4 ) 

/CSEVAL/, 
17, 22-25, 
35, 37-41, 
77 
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TABLE 5. SUBROUTINE BARSET NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

1)ELT 

Temperature! difference 
between two values in C -T 
tabic! 

0 R 

24, 25, 40, 
70, 82 

ERASE 1 

Minus value of stagnation 
specific! he!at multiplied by 
the! natural logarithm of 
TCTAB ( !) 

Btu/(lbm-° It) 

71, 72, 80 

FJ 

Conversion factor between 
thermal and work units 

( ft-lbf) /Btu 

/INPUT/, 
2, 2, 59 

FJG 

FJ multiplied by G 

(ftj-lbm) 

(Btu-sec 2 ) 

/CSEVAL/, 
2, 60 

G 

Acceleration of gravity 
used as a proportionality 
constant 

(lbm-ft) 

(lbf-sec 2 ) 

TFnput/, 

2 

G1 

Summation of entropy 
terms 

Btu/ ( lbm-° R) 

26, 40, 42 

G2 

Intermediate term for 
entropy equation 

Btu/( lbm-° R) 

41, 42 

GAMO 

Specific heat ratio at 
stagnation condition 


/INPUT/, “ 

46-48, 50, 

51 

GM102 

One-half the specific 
heat ratio for stagnation 
condition minus one 


/CSEVAL/, 

50 

GOG Ml 

Specific heat ratio at 
stagnation condition 
divided by the specific heat 
ratio minus one 


/CSEVAL/, 
51, 59 

GTAB 

Array of terms used by 
SEVAL 

Btu/( Jbm-° R) 

/CSEVAL/, 
25, 42, 73 

HO - 

Stagnation enthalpy 

ft*/sec^ 

/CSEVAL/, 
45, 61 

HTAB 

Array of enthalpy values 

Btu/lbm 

/CSEVAL/, 
21, 35, 81 

I 

Do loop counter 


8-10, 27, 31, 
35-42, 53, 

54, 62, 63, 
72-83 
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TABLE 5. SUBROUTINE BARSET NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

II 

Saved value of one or 
IXTAB 

— 

5, 7, 8 

ICTAB 

Indicator; ~ 0, constant 
specific heat calculation; 
= (3 ICTAB ^ 20), 

number of C values in 
P 

input table 


/INPUT/, 
12-14, 16, 
17, 27, 36, 
44, 55 

IS 

Number of C values in 
P 

input table minus one 

mB 

/CSEVAL/, 
44, 58 

ITWTAB 

Indicator; = -1, T 

w 

* adiabatic wall tempera- 
ture; = 1, IXTAB input 

values of T are used 
w 


/INPUT/, ~ 

6, 85 

EXTAB 

Number of points in the -x, 
y, and Mach number tables 

— 

/INPUT/, 
7, 53, 62 

NOCTAB 

Saved value of ICTAB 

1 

/CSEVAL/, 
16, 56, 57, 
67, 70, 72 

NOCTM1 

NOCTAB minus one 

— 

1 

PO 

Stagnation pressure 

mm 

/INPUT/, 
52, 87 

POMAX 

Saved value of stagnation 
pressure 

m/it 2 

/CSEVAL/, 

52 

PIE 

Circular constant tt 

— 


PITAB 

Pressure table for nozzle 
contour points, x, y, 
corresponding to IXTAB 
values; obtained from 
GETPT 

lbf/ft 2 

/TABLES/, 

54 

RBAR 

Specific gas constant: 
universal gas constant 
divided by molecular 
weight 

■ 

/INPUT/, 
3, 59 
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TABLE 5. SUBROUTINE BARSET NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

ROJ 

Specific gas constant in 
thermal uni Is , Hi iA J v 
divided by FJ 

Btu/( lbm-° R ) 

/CSEVAL/, 
3, 46 

SO 

Stagnation entropy 

Btu/ ( lbm- c It) 

/CSEVAL/, 

87 

TO 

Stagnation ( empcralure 

° R 

/INPUT/, 
4, 45, 61, 
87 

TCTAB 

Temperature values 

corresponding to ICTAB 

values in increasing order 

defining C -T tables 
P 

°R 

/CSEVAL/, 
13, 14, 17- 
20, 31, 67, 
68, 70, 71, 
81, 83 

TE 

Saved value of TITAB 

°R 

63-65 

TIE1 

Saved value of TCTAB (I) 

°R 

18, 21-25, 
28, 31-34, 
37-41 

TIE 2 

Squared value of TIE1 

“R ! 

19, 21-23, 
29, 33, 37, 
38, 40, 41 

TIE 3 

Cubed value of TIE1 

or 3 

20, 22, 30, 
33, 37, 40, 
41 

TITAB 

Temperature table 
corresponding to IXTAB 
values for contour points x, 
y; obtained from GETPT 

°R 

/TABLES/, 
54, 63 

TMAX 

Saved value of maximum 
temperature 

° it 

4, 9, 10, 
13, 14, 64, 
65, 67 

TME1 

Saved value of TIE1 

°R , 

28, 34, 40 

TME2 

Saved value of TIE 2 

°R 2 j 

29, 40 

TME3 

Saved value of TIE 3 

°R S 

30, 40 

TWTAB 

Wall temperature table 
corresponding to IXTAB 
values for wall locations 
x, y. 

°R 

/TABLES/, 
9, 10, 86 

Z1 

Real value of NOCTM1 


69, 70 



SUBROUTINE BMFITS 


Subroutine BMFITS computes the coefficients from a cubic spline 
fit routine using table la r data. 

COMMON BLOCKS . . 

No COMMON blocks arc used. 

TBL SUBROUTINES 

Subroutine BARSET calls BMFITS. 

BMFITS does not call any TBL subroutines. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines or built-in FORTRAN functions are 
used. 

CALLING SEQUENCE 

The subroutine calling sequence is 
CALL BMFITS (X, Y,N, BL,CL,DL) 

where 

X = array of independent table values, 

Y ~ array of dependent table values, 

N = number of values in the X and Y arrays, 

BL = array of coefficients for first-power term, 

CL = array of coefficients for second-power term, 

DL = array of coefficients for third-power term. 
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SOLUTION METHOD 


Initialize array subscript counter to one. 

1. I = 1 

Compute difference between adjacent values in X array. 

2. FL(I) = X(I + 1) - X(I) 

Increment subscript counter. 

3. 1=1+1 

Check whether subscript counter has covered the complete X array, 

4. If I < N, go to 2 
If I a N, go to 5 

Set subscript counter to two. 

5. 1=2 

Compute intermediate terms. 

6. B(I) = ~(FL(I - 1))/(FL(I)) 

7. A(I) = -(2.0)(FL(I) + FL(I - l))/(FL(l)) 

8. C(l) = ((6.0)/(FL(l))((Y(l + 1) - Y(I))/(FL(I)) - (Y(l) 

-Y(I - 1))/(FL(I - 1))) 

Increment subscript counter. 

9. 1=1+1 

Check whether subscript counter has covered the complete array. 

10. If I < N, go to 6 
If I ^ N, go to 11 



Initialize intermediate terms. 


11. G(2) = 1.0 

12. F(2) = 0.0 

Set subscript counter to three. 

13. 1=3 

Compute intermediate terms. 

14. G(I) = A(I - 1) + (B(I - 1))/(G(I .. 1)) 

15. F(I) = -( (B(l - 1))(F(I - 1))/(G(I - 1)) + C(l - 1)) 
Increment the subscript counter. 

16. 1=1+1 

Check whether the subscript counter has covered the complete array. 

17. If I ^ N, go to 14 
If I > N, go to 18 

Compute intermediate array terms starting with the last value. 

18. YPP(N) = (F(N))/(G(N) - 1.0) 

19. YPP(N - 1) = YPP(N) 

20. I = N - 2 

21. YPP(l) = (YPP(I + 1) + F ( I + 1))/(G(I + 1)) 

Decrement subscript counter. 

22 . 1 = 1-1 


Check whether subscript counter has reached zero. 


23. If I > 0, go to 21 
If I — 0, go to 24 

Set subscript counter to one. 

24. 1=1 

Compute array of polynomial coefficients. 

25. BL(I) = (Y (i + 1) _ Y(i))/(FL(I)) - ( (FL(l) ) (YPP(l + l) 

+ (2.0)(YPP(l))))/(6.0) 

26. CL(I) = (YPP(l))/(2.0) 

27. DL(l) = (YPP (i +1) - YPP(I))/((6.0)(FL(I))) 

Increment subscript counter. 

28. 1=1+1 

29. If I < N, go to 25 
If I — N, go to 30 

30. Return 

Table 6 gives subroutine BMFITS nomenclature. 


TABLE 6. SUBROUTINE BMFITS NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

A 

Array of intermediate terms 

— ■ 

DIM, 7, 14 

B 

Array of intermediate terms 

— 

DIM, 6, 14, 
15 

BL 

Output array of polynomial 
coefficients 

Btu/(lbm-°R 2 ) 

CALL, DIM, 
25 

C 

Array of intermediate terms 

— 

DIM, 8, 15 

CL 

Output array of polynomial 
coefficients 

Btu/( lbm-° R 3 ) 

CALL, DIM, 
26 
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TABLE 6. SUBROUTINE BMFITS NOMENCLATURE (Concluded) 


Symbol 


Description 


Units 



DT. Output array of polynomial Btu/( lbm-° R 4 ) 

coefficients 

F Array of intermediate terms — 

FL Array of intermediate terms 0 R — — 

G Array of -intermediate terms - — 


Array subscript counter 


Input dimension of table 
arrays 


Input array of independent 
table values 


Input array of dependent 
table values 


Array of intermediate terms 



Reference 

CALL, DIM, 

27 

'DIM, 12, 

15, 18, 21 
DIM, 2, 6-8, 

25, 27 

DIM, 11, 

14, 15, 18, 

21 


1-10, 13-17, 
20-29 


CALL, 4, 
10, 17-20, 
29 


CALL, DIM, 
2 


CALL, DIM, 
8, 25 


DIM, 18, 19, 
21, 25-27 














SUBROUTINE DIRECT 


Subroutine DIRECT controls the overall flow of the program. It calls 
subroutines which successively read in data, compute program constants, fit 
curves, and then obtain the boundary layer solution. 

COMMON BLOCKS 

No COMMON blocks arc used. 

TBL SUBROUTINES 

Subroutines MAINTB and QUITS call DIRECT. 

DIRECT calls BARCON, BARSET, and READIN. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines or built-in FORTRAN functions are 

used. 

CALLING SEQUFNCE 

The subroutine calling sequence is: CALL DIRECT. 

Call subroutine to read input data. 

* 

1. CALL READIN 

Call subroutine to compute constants and initialize data. 

2. CALL BARSET 

Obtain the boundary layer solution . 

3. CALL BARCON 
Read in the next case. 

4. Go to 1 

t 

! 

i 

% 

i 
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SUBROUTINE CFEVAL 


Subroutine CFEVAL evaluates the skin friction coefficient cT f as a 
function of C f Rj based on Coles' relationship. For values of C f R- such 
that 0 < R~ < 2,01, C is evaluated from 

,, _ 0. 009890 

f (C f l^) 6 * 882 • 

For values of C R 3 , such that 2.51 =£ C R < 1982759.2, C is 

1 0 I 6 f 

evaluated from a cubic spline fit of log (C ) versus log (C. R^-) from input 

1 10 

data. All other values of C^ R^ are considered out of range. 

COMMON BLOCKS 

No COMMON blocks are used. 

TBL SUBROUTINES 

Subroutines BARPRO and START call CFEVAL. 

CFEVAL calls subroutine QUITS, 

FORTRAN SYSTEM ROUTINES 

FORTRAN library routines A LOG and EXP are used. 

No built-in FORTRAN functions are used 
CALLING SEQUENCE 

CFEVAL is a function routine and is called as: 

CFBAR = CFEVAL (CFRT), 

where CFRT is the skin friction coefficient multiplied by the Reynolds 
number based on the momentum thickness, and CFBAR is the skin 
friction coefficient. 

SOLUTION METHOD 

Save r R fl - = (CFRT) . 

I V 

1. Z = CFRT 
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Check the sign of C R_, 
t 6 


2. If Z JS 0, go to 6 
If Z > 0, go to 3 

Check whether the integer IZ of Z * (C^-) greater, equal, or 

smaller than the integer IX(J) of the X table values. 

3. If IZ > IX(J), go to 18 
If IZ = IX( J), go to 13 
If IZ < IX( J), go to 4 

Set 

4. J = J - 1 

Check 

5. If J > 0, go to 3 
If J = 0, go to 9 
If J < 0, go to 6 

Set J equals one and print a message. 

6. J = 1 

7^- -WRITE Z, ZERO, X(8) 

Stop calculation and go to next case. 

8. CALL QUITS 

Check the sign of Z = CFRT. 

9. If Z s 0.0, go to 6 
If Z > 0.0, go to 10 


Set 

10. J = 1 

„ , ~ y,,v , .. , 77 0.009896 

Compute (Y) from the relation = (g 1 ^ yi.siia 

11. Y = (0. 009896) /(Z)°* 582 

12. Go to 16 

Save the natural logarithm of Z = CFRT, 

13. ZL = Ln(Z) 

Compute ln(C^) = YL by using the polynomial equation obtained from 
Coles’ experimental data, in the case 2.51 s C , Kr < 1982759.2. 

I V 

14. YL. = D(J) + (ZL)(C(J) + (ZL)(B(J) + (ZL)(A(j)))) 

Obtain C^ = Y from -ln(C t .) . 

15. Y = Exp (YL) 

Define CFEVAL representing C^. 

16. CFEVAL = Y 

17. Return 

Compare the integers of IZ and IX (J + l). 

18. If IZ £ IX( J + 1), go to 13 
If IZ > IX( J + 1) , go to 19 

Set 

19. J = J + 1 

Check whether J exceeds eight. 
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20. If J a 8, go to 6 



If J < 8, go to 18 

Table 7 gives subroutine CFEVAL nomenclature. 


TABLE 7. SUBROUTINE CFEVAL NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

A 

Coefficient in ln(cp versus 

— 

DATA, DIM, 14 


In (Oj, K-) polynomial 



B 

Coefficient in In fC^) versus 

— 

DATA, DIM, 14 


In (C„ R s ) polynomial 
t o 



C 

Coefficient in In (C f ) versus 

— 

DATA, DIM, 14 


In (C f R^) polynomial 



CFEVAL 

Saved value of Y, skin friction 
coefficient C 

— 

16 

CFRT 


— 

CALL, 1 

D 

Coefficient in In (U^) versus 
In (C. R-) polynomial 

I v 


DATA, DIM, 14 

IX 

i 

Integer value of X 



DIM, EQUIV, 3, 
18 

IZ 

Integer value of Z 

— 

EQUIV, 3, 18 

.T 

Indicator to pick up data 

— 

DATA, 3-6, 10, 
14, 18-20 

X 

Input value of C, H s (data) 

I C7 

— 

DATA, DIM, 
EQUIV, 7 

Y 

C { , Skin friction coefficient 

— 

11, 15, 16 

YL 

In (C { ) 

— 

14, 15 

Z 

Saved value of CFRT = C^ R- 

— 

EQUIV, 1, 2, 7, 
9, 11, 13 

ZERO 

0.0 

_ 

DATA. 7 

ZL 

to <e, %> 

— 

13, 14 
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SUBROUTINE FIIF 


Subroutine FIIF defines the function (of S) to be evaluated by the 
numerical integration routine INTZET. Three second-degree coefficients 
(Ap, B F , C p ) and an exponent value n are determined by subroutine ZETAIT. 

Enthalpy value h is computed from 


h 


= H + 
w 



S + 



= A, 


b p s 


c F si 


(AFINT) (BFINT) (CFINT) 


Entering subroutine SEVAL with h, the temperature t is obtained. 
The appropriate form of the function is then evaluated for 


/t\ 

IFINT = 1, f (s) = yj* I 


S* 1 (1 - S) 


• «i 


FIIF (S) 


and 


IFINT = 2, f (s) =( Y J s ••• FIIF (S) 


where T /F has been used for p/p 
e ^ r e 


COMMON BLOCKS 

COMMON block COFUE_is used. 
TBL SUBROUTINES 


Subroutine INTZET calls FIIF. 


FIIF calls SEVAL. 


FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines or built-in FORTRAN functions are 
used. 


S8 


CALLING SEQUENCE 

FIIF is a function routine and is called as 
Y = FIIF(S) , 

where S is an independent variable 0 s S — 1 and defined as 



FIIF (s) represents 

fs”(l-S).c^ . 

SOLUTION METHOD 

Initialize power summation term to 1.0. 

1. STOM = 1.0 

Check the sign of power index n used in the velocity and enthalpy 
distribution equations. MMINT = (n) is defined in subroutine ZETAIT. 

2. If MMINT ^ 0, go to 5 
If MMINT > 0, go to 3 

Save ,! S n " as STOM. 

3. Do 4. I = 1, MMINT 
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4. STOM = (STOM) (s) 
Calculate the enthalpy. 


h = H + 
w 


H - H 
o w 


( 4 ) 


s 2 . 


5. FDEN = AFINT + (BFINT) (S) + (CFINT) (s ) 2 
Check the option indicator, 
t. If IFINT 2= 2, go to 9 
If IFINT < 2, go to 7 
Calculate S n (l - S). 

7. FNUM = (STOM) (1.0 - S) 


8. Go to 10 


_ , n 

Save & . 


9. FNUM = STOM 

Call subroutine SEVAL to obtain the temperature t = (T) and specific 
heat C p = (<j>) by using the previously determined enthalpy 

h = (FDEN). 

10. CALL SEVAL (2, T, O, FDEN) 

Calculate the function FIIF. 

11. FIIF = (FNUM)(TFINT)/(T) 

Where TFINT is the free-stream temperature T determined in 
BARPRO and defined in ZETAIT, e 

12. Return 

Table 8 gives subroutine FIIF nomenclature. 
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TABLE 8. SUBROUTINE FIIF NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

AFINT 

Coefficient in enthalpy 
equation representing 
H w ( wa ^ enthalpy) . 

ft 2 /sec 2 

/ COFIIF/, 5 

BFINT 

Coefficient in enthalpy 

equation representing 

(H - H )/£, 
v o w 

ft 2 /sec 2 

/ CO FIIF/, 5 

CFINT 

Coefficient in enthalpy 
equation representing 

(-UV2). 

ft 2 /sec 2 

/COFIIF/, 5 

FDEN 

Enthalpy 

h = H + E (H - H )/£]S 
w o w 

- (U 2 /2)S 2 
e 

ft 2 /sec 2 

5, 10 

FNUM 

n „n. 

S or S (1 - S) 

, 

7, 9, 11 

IFINT 

= 1; Calculation of 

fs n d - s) 

= 2; Calculation of 
T 

-^s n 

t 


/COFIIF/, 6 

MMINT 

Velocity power law 
exponent n 

— 

/COF-HF/, 2, 
3 

~0 

Specific heat C 

P 

Btu/( lbm-° R) 

10 

S 

- / \ 1//n 
■■f w * 

h - H / \l/n 

e_o w/y\ 

H - 11 \A / 

0 w ' ' 


CALL, 4, 5, 
7 


101 















SUBROUTINE GETPT 


Subroutine GETPT computes pressure P c 
H, and T relations for a given Mach number 


and temperature T^ from 

M . 

e 


For constant specific heat y^ = y Q 



For variable C the following iteration is used. 

P 

1. An approximation is made for temperature T 

eg 


using 


T = 
eg 


- 1 


1 + 


M‘ 


2. In subroutine SEVAL, 


C and H are computed at temperature 
pe e 


eg* 

C = A + BT + CT 2 + DT 3 ; 
P 


H(T) = H. + / C (t)dt . 
1 T. P 


3. The specific heat ratio y is then evaluated from 


C 



4. Tempei'aturc T is then calculated, 
c 


T 

e 


2(H - H ) 

o e 


y RM 
e e 


since 

a 2 = y RT , 
e e 


U 2 U 2 

y RM 2 = y R -J = — 1 , 

q o ' o 06 nn ? 


U 2 

T e 

e y RM 2 * 


U 2 

H = H + — - , 

o e 2 


U 2 = 2(H - H ) = T y RM 2 

e o e e e e 


T = 


2(H - H ) 


e ~y RM 2 
e e 


5. A relative error comparison is then made. If 


T - T 
e eg 

T 

eg 


^ TOLZME 


9 


convergence is obtained, and subroutine SEVAL is used to 

evaluate pressure P at temperature T . If the convergence 

e e 

criterion is not satisfied, a form of Wegstein* s method is used to 
approximate T again, and steps 2 through 5 are repeated. A 

maximum of 50 iterations are possible. 


COMMON BLOCKS 


COMMON blocks CSEVAL and. INPUT aroused. 

TBL SUBROUTINES 

Subroutines BARPRO, BARSET, and START call GETPT. 
GETPT calls subroutine SEVAL. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines are used. 

Built-in FORTRAN function ABS is used. 

CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL GETPT (ZME, PI, TI) 

where 

ZME = free -stream Mach number, 

PI = free-stream pressure, 

TI = free-stream temperature. 

SOLUTION METHOD 

Determine the square of Mach number. 

1. ZME2 = (ZME) 2 

Calculate 

gRM 2 . 


2. PROD1 = (2.)/(RBAR)/(ZME2)/(G) 


Compute: 



2 



3. DENM2 -- 1.0+ (GM102)(ZME2) 


Obtain froe-stream temperature 



4. TK - (T0)/(DENM2) 

Check the value of ICTAB [= 0: constant specific heat calculation, 
= (3 ^ ICTAB =£ 20): number of points in variable C versus T 
table] . 

5. If ICTAB > 0, go to 10 
If ICTAB =£ 0, go to 6 

Calculate 


P 

e 




y /y -l 

e e 


for y 

e 


constant. 


6. PE = (P0)/(DENM2) GOGM1 


Set 


7. PI = PE 


8. TI = TE 

Return to the main routine or continue the calculation unless 

C is constant, 
pe 

9. Return 





Set the iteration counter to zero initially. 


10. ITER = 0 
Define the tolerance. 

11. TOL ~ (TOLZME) / (ZME) 

Save the previous approximation of free -stream temperature T . 

eg 

12. TEO = TEG 

2 (II - H ) 

O G 

Save T^ = — — as calculated in step 17. 

c ^e 1 e 

13. TCO = TC 

Approximate T . 

eg 

14. TEG = TE 

Call subroutine SEVAL to obtain C and H for the known T . 

pe e e 

15. CALL SEVAL(l, TE, CPE , HE ) 

Calculate specific heat ratio y g . 

16. GAME = (CPE) /(CPE - ROJ) 

2 (H - H ) 

0 e 

Calculate new free-stream temperature • 

c *e e 

17— TC = (HO - HE)(PRODl)/(GAME) 

Test whether the temperature falls within the tolerance. 

18. If I (TC - TE)/(TE)| ^ TOL, go to 33 
If I (TC - TE)/(TE) I > TOL, go to 19 
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Check on the number of iterations. 


19. If ITER > 0, go to 22 
If ITER 0, go to 20 

Determine 

20. TK = ((2.0) (TK) + TC)/(3.0) 

21. Go to 30 

Check on the number of iterations, and print out the error message 
if convergence was not obtained within 50 iterations, 

22. If ITER ^ 50, go to 25 
If ITER > 50, go to 23 

Print out error message, 

23. WRITE ZME, TC, TCO, TE, TEO 

24. Go to 33 

Use a form of Western's method to approximate a new value for T . 

e 

25. ZK = (TC - TCO)/(TE - TEO) 

26. TE = (TC - (ZK)(TE))/(1.0 - ZK) 

Check whether convergence has been achieved. 

27. "If I (TE - TEG)/(TE) | < TOL, go to 32 

If | (TE - TEG)/(TE) | * TOL, go to 28 
Check on the number of iterations, 

28. If ITER < 10, go to 30 
If ITER 2= 10, go to 29 

Check whether the calculated number falls within the tolerance. 
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29. If I (TK - TEO)/(TE)| < TOL, go to 92 
If | (TK - TKO) / (TK) | & TOL, go to 90 
Add ono to tho iteration counter. 


DO. ITER ITER + 1 

Repeat tilt* abovo method until a convergence Is obtained. 

31. Go to 12 

Determine the new free-stream temperature. 

32. TK « (TE + TEG) / (2. 0) 


Call subroutine SEVAL to obtain the free-stream pressure P^, 


using temperature T^ and stagnation entropy S^. 

33. CALL SFVAL (-1, TE, PE, SO) 


34. Go to 7 

Table 9 gives subroutine GETPT nomenclature. 

TABLE 9. SUBROUTINE GETPT NOMENCLATURE 


Symbol 


CPE 


DENM2 


Description 


Specific heat at con. 
stant pressure in 
free-stream C 

pe 


7-1 

1 + M 2 

2 e 


Acceleration of 
gravity used as a 
proportionality 
constant 


Units 


Btu/( lbm- e R) 


( lbm/lbf ) /( ft/ sec 2 ) 


Reference 


15, 16 


3, 4, 6 


/INPUT/, 2 


* 

t 


















1 ABLE 9, SUBROUTINE GETPT NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

GAME 

Specific heat ratio 
in fret! stream 

— 

1(5, 17 

GM102 

2 

— 

/CSEVAL/, 3 

GOGM1 

y e 

Y e - 1 


/CSEVAL/, 6 

HO 

Stagnation enthalpy 
H 

o 

ft 2 /sec 2 

/CSEVAL/, 

17 

HE 

Free -stream 
enthalpy H 

e 

ft 2 /sec 2 

15, 17 

ICTAB 

Indicator; - 0, con- 
stant specific heat 
calculation; 

= (3 S ICTAB ^ 20), 

dimension of variable 

C versus T table 
P 


/INPUT/, 5 

ITER 

Iteration counter 


10, 19, 22, 
28, 30 

RO 

Stagnation pressure 
P 

o 

lbf/ft 2 

/INPUT/, 6 

PE 

Free -stream 
pressure P 

e 

lbf/ft* 

6, 7, 33 

PI 

Saved value of PE 

IbfTft 2 

CALL, 7 

PROD1 

2 

gRM^ 

(°R-sec 2 )/ft 2 

2, 17 

RBAR 

Specific gas con- 
stant; universal 
gas constant 
divided by molec- 
ular weight 

( ft-lb) / ( lbm-° R) 

/INPUT/, 2 


























TABLE 9. SUBROUTINE GETPT NOMENCLATURE (Concluded) 


Description 


Reference 


TOLZME 


ZME 


ZME2 


Specific gas constant Btu/(lbm-*R) 
in thermal units, 

RBAR divide d by FJ 

Stagnation entropy Btu/( lbm-° R) 








Stagnation tem- 
perature T 

o 


Saved value of 
free -stream 
temperature 


Saved value of TC 


Free-stream 
temperature T 


Approximation of 
free-stream tem- 
perature in iteration 
loop 


Saved value of TEG 


Saved value of TE 


Tolerance value 
divided by free- 
stream Mach number 


Tolerance 


Intermediate term 
in free-stream tem- 
perature equation 


Free-stream Mach 
number M 



/CSEVAL/, 

16 

/cseval/, 

33 

/INPUT/, 4 


13, 17, 18, 

20, 23, 25, 26 


13, 23, 25 


4, 8, 14, 15, 
18, 20, 23, 
25, 26, 27, 
29, 32, 33 


12, 14, 27, 
32 


12, 23, 25, 

29 


CALL, 8 


11, 18, 27, 
29 


/INPUT/, 11 


25, 26 


CALL, 1, 11, 23 


1, 2, 3 



I 






















SUBROUTINE INTZET 



%'■ 


Subroutine INTZET uses a numerical technique in connection with a 
quintic spline fit to evaluate an integral within the limits X t = XC(l) and 
X 2 = XC(2l) . The integration interval is divided into 20 coarse subintervals. 
Since the integration limits are restricted to certain values (0, 1, £, l/£ ), a 
test on the interval distance (X 2 - Xj) is performed. If (X 2 - Xj) =s o. 1, the 
function is evaluated at 21 subinterval endpoints, and the integral is approxi- 
mated by the quintic technique with Simpson’s rule for the outer subintervals. 
Eor (X 2 - X*) > 0. 1 the maximum function value within the interval is deter- 
mined, and an approximation of the integral is performed over the coarse 
subintervals up to a point where the function value reaches 80 percent of the 
maximum value. Beyond this point the remaining coarse subintervals are 
further subdivided individually into three fine subintervals. Solving the function 
at the endpoints of the fine subintervals, the integral approximation is finally 
completed using these finer subinterval step sizes. 

For a given function f(s), evaluated at equally spaced (As) points 
S o , Sj, S 2 , s 3 , the quintic approximation of the integral is prepresented by 

r As 

f f(S)dS = — [-f(S Q ) + 13f(S t ) + 13f(S2) - f(S 3 )] . 

The one-sided Simpson rule is expressed by the following equation 
and diagram. 

1 

/ f(s)dS = [5f(S Q ) + 8f(S t ) - f(S 2 ) ] . 
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COMMON BLOCKS 


No COMMON blocks are used. 

TBL SUBROUTINES 

Subroutines START and ZETAIT call INTZET, 
INTZET calls FIIF. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines are used. 
Built-in FORTRAN function FLOAT is used. 
CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL INTZET (XI, X2, ZINT ) 
where 

XI = lower limit of integration 
X2 = upper limit of integration 
ZINT = integral value. 

SOLUTION METHOD 

Set the length of interval for integral. 

1. DX21 = X2 - XI 

Set the initial value of integral as zero. 

2. SUMINT =0.0 

Check whether the interval length is zero. 

3. If DX21 = 0.0, go to 22 
If DX21 * 0.0, go to 4 


Divide the interval over which the integration is performed into 20 
coarse subintervals. 

4. DXC = ( DX2l) /( 20. 0) 

Set maximum value of integral and its subscript counter to large 
negative values. 

5. I MAX = -9999 

6. FMAX - -1.0E30 

Determine the independent variable, the integrand of the function, and 
the maximum value of the integrand. 

7. Do 14, I = 1, 21 
Independent variables 

8. XC(I) = XI + (FLOAT(l - l))(DXC) 

Call subroutine FIIF to determine the integrand. 

9. YC (i) = FEF(XC(I)) 

Check for the maximum value of the integrand. 

10. If YC(l) ^ FMAX, go to 14 
If YC(l) > FMAX, go to 11 

Save subscript counter at the maximum value of the integrand. 

11. IMAX = I 

Save independent variable corresponding to maximum integrand. 

12. XMAX = XC(I) 

Save maximum value of integrand. 

13. FMAX = YC (i) 
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14. Continue 


Check for the interval over which the coarse or fine integration must 
be performed. 

15. If DX21 > 0.10, go to 24 
If I)X21 ^ 0. 10, go to 10 

Determine 24 times the quantity of the integral according to the one- 
sided Simpson' s Rules 

16. SUMINT = ( 10. 0) ( YC( l) ) + ( 16. 0) ( YC( 2) ) - (2.0)(YC(3)) 
Calculate the quintic approximation of Simpson' s rule. 

17. Do 19, I = 2, 19 

18. PARINT = ( 13. 0) ( YC(l) + YC(l + l)) - YC(l - l) - YC(l + 2) 

19. SUMINT = SUMINT + PARINT 

20. SUMINT = SUMINT + (10.0) (YC(2l) ) + (16.0) (YC(20))- 

(2.0)(YC(19)) 

Obtain the integral. 

21. SUMINT = (SUMINT)/(24.0)(DXC) 

Save the integral 

22. ZINT = SUMINT 
Go back to main routine. 

23. Return 

Determine the function value 20 percent off the maximum value. 

24. FBRK = (FMAX) (0. 20) 
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Set the initial value of the integral equal to zero, 

25. SUAINT =0.0 

26. SUBINT =0.0 

Check the value of the subscript counter for the maximum value of the 
integrand. 

27. If IMAX s 2, go to 35 
If IMAX > 2, go to 28 

Determination of IBRK corresponding to FBRK is made from step 28 
to step 33. 

28. Do 32, I = 2, IMAX 

Check whether the integrand has reached the switch-over value. 

29. If YC(l) ^ FBRK, go to 32 
If YC(l) > FBRK, go to 30 

Save the subscript counter for the switch-over value. 

30. IBRK =1-1. 

31. Go to 34 

32. Continue 

Set the swich-over subscript counter to one less than the maximum. 

33. IBRK = IMAX - 1 
Check the value of IBRK. 

34. If IBRK > 1, go to 38 
If IBRK ^ 1, goto 35 
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Set the switch-over subscript counter and switch-over indicator, 

35. IBRK = 1 

36. IBRKM1 = 0 

37. Go to 45 

Save IBRK minus one. 

38. IBRK Ml = IBRK - 1 

Calculate 24 times the quantity of the integral according to the one- 
sided Simpson* s rule. 

39. SUAINT = (10.0)(YC(1)) + ( 16. 0) ( YC( 2) ) - (2.0)(YC(3)) 
Check the value of IBRK. 

40. If IBRK ^ 2, go to 44 
If IBRK > 2, go to 41 

Continue integration. 

41. Do 43 I = 2, IBRKM1 

Compute 24 times the quantity of the integral for each interval using 
Simpson* s rule. 

42. PARINT = ( 13. 0) ( YC(l) + YC(l + l)) - YC(l - l) - YC(l + 2) 
Sum the individual interval integrals. 

43. SUAINT = SUAINT + PARINT 

Obtain the integral up to where the integrand is 20 percent of the 
maximum value. 

44. SUAINT = (SUAINT) (DXC)/( 24.0) 

Obtain one-third DXC to consider finer subintervals. 
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45. UXM = (DXC) /(;}, 0) 

Check the value of IBRKM1. 

46. If IBRKM1 > 0, go to 50 
If IBRKM1 S 0, go to 47 

Set temporary integers. 

47. K = 2 

48. JS = 2 

49. Go to 52 

Set temporary integers. 

50. K = 3 

51. JS = 1 

Determine the independent variables corresponding to fine subintervals. 

52. Do 54, I - 2, 4 

53. XM = XC(IBRK) + (FLOAT(l - K))(DXM) 

Call subroutine FIIF to obtain integrand. 

54. YM(I) » FIIF(XM) 

Check the value of IBRKM1. 

55. . If IBRKM1 > 0, go to 57 

If IBRKM1 ^ 0, go to 56 

Calculate 24 times the quantity of the integral corresponding to the fine 
subinterval. 


56. SUBINT = (10.0)(YM(2)) - (2.0)(YM(4)) + ( 16 . 0) (YM(3)) 



Continue the integration dividing the course subintervals into three 
fine subintervals. 

57. Do 07, I - TURK, If) 

58. Do 05, eJ .IS, ;> 

59. YM( i) - YM(2) 

00. YM(2) - YM(:t) 

61. YM(.*i) = YM(4) 

62. XM = XM + DXM 

Call subroutine FIIF to obtain integrand. 

63. YM(4) = FIIF(XM) 

Compute 24 times the. quantity of the integral for fine subintervals. 

64. PARINT = (13.0)(YM(2) + YM(3)) - YM(l) - YM(4) 

Sum the subinterval integrals. 

65. SUBINT = SUBINT + PARINT 
Set 

66. JS = 1 
Independent variable; 

67. XM = XC(I + 1) + DXM 

So far, 24 times the quantity of the integral up to XC(20) has been 
obtained. Calculate the remaining integral between XC(20) and 
XC(2l). 

68. Do 75, J = 1, 2 

69. YM( 1) = YM(2) 
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70. YM(2) = YM(3) 

71. YM(n) = YM(4) 

72. XM - XM + DXM 

Call .subroutine FIIF to obtain integrand. 

72. YM(4) = FIIF(XM) 

74. l’A HINT = (13.0)(YM(2) + YM(a)) - YM(l) - YM(4) 

Calculate 24 times the quantity of the integral up to the independent 
variable of XC(20) + ( 2.0) (DXM). 

75. SUBINT - SUBINT + PARINT 

Determine 24 times the quantity of the integral up to XC(2l). 

76. SUBINT = SUBINT + (10. 0) (YM(4) ) + (l6.0)(YM(3) 

- (2.0)(YM(2)) 

Calculate the final integral beyond the point where the value of the 
integrand is 20 percent of the maximum value. 

77. SUBINT = (SUBINT) (DXM)/(24.0) 

Obtain the total integral. 

78. SUMINT = SUAINT + SUBINT 

79. Go to 22 

Table 10 gives subroutine INTZET nomenclature. 


120 


sr 


TABLE 10. SUBROUTINE INTZET NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

DX21 

Upper limit of integration 
minus lower limit 


1, 3, 4, 15 

DXC 

Coarse subinterval 

— 

4, 8, 21, 44, 45 

DXM 

Fine subinterval 

— 

45, 53, 02, 07, 
72. 77 

FBRK 

20% value of maximum 
integrand 

— 

24, 29 

FMAX 

Saved value of maximum 
integrand 

— 

6, 10, 13, 24 

I 

Do-loop counter 


7-13, 17, 18, 
28-30, 41, 42, 
52-54, 57 


Subscript counter for 
switch-over value of the 
integrand 


30, 33-35, 38, 
40, 53, 57 

IBRKMl 

IBRK- minus one 

— 

36, 38, 41, 46, 
55 

IMAX 

Subscript counter value 
corresponding to maxi- 
mum integrand 


5, 11, 27, 28, 
33 

J 

Do -loop counter 

— 

58, 68 

JS 

Starting do -loop value 

— 

48, 51, 58, 66 

K 

Switch-over point 
counter 

■■■■ 

47, 50, 53 

PARlNT 

24 times the quantity of 
integral 

— 

18, 19, 42, 43, 
64. 65. 74. 75 

SUAINT 

Integral value up to 20% 
of maximum integrand 

— 

25, 39, 43, 44, 
78 


Integral value beyond 
20% of maximum 
integrand 


26, 56, 65, 75-78 

SUMINT 

Integral 

— 

2, 16, 19-22, 78 

Xi 

Lower limit of 
integration 

— 

CALL, 1, 8 

X2 

Upper limit of 
integration 


CALL, 1 

XC 

Independent variable 


DIM, 8, 9, 12, 
53, 67 


1 t 
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TABLE 10. SUBROUTINE INTZET NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

XM 

Independent variable 
corresponding to fine 
subintervals 

— 

53, 54, 82, 63, 
67, 72, 73 

XMAX 

Saved value of the 
independent variable 
corresponding to maxi- 
mum value of integrand 


12 

Y~C 

Integrand 


DIM, 9, 10, 13, 
16, 18, 20, 29, 
39, 42 

YM 

Integrand corresponding 
to fine subintervals 


DIM, 54, 56, 59- 
61, 63, 64, 69- 
74, 76 

ZINT 

Saved value of integral 

— 

CALL, 22 
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SUBROUTINE MAINTB 

Subroutine MAINTB provides the main control for the turbulent 
boundary layer program. 

COMMON BLOCKS 

COMMON block INPUT is used. 

TBL SUBROUTINES 

MAINTB calls subroutine DIRECT. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines or built-in FORTRAN functions 
are used. 

SOLUTION METHOD 

Initialize step -size indicator. 

1. IDXMAX = 0 

Call the overall flow control program. 

2. CALL DIRECT 

Table 11 gives the nomenclature for subroutine MAINTB, 


TABLE 11. SUBROUTINE MAINTB NOMENCLATURE 


Symbol 


Description 


Units 


Reference 


IDXMAX 


Indicator for previous value 
of step size or to compute the 
step size 

= 0, compute step size; 

= 1, use previous step size 


/INPUT/, 1 


1 

I 

i 

1 






SUBROUTINE QUITS 


Subroutine QUITS is entered when an error has been detected in the 
input or if an unreasonable number has been calculated during execution. It 
prints out appropriate error statements and the contents of the COMMON 
blocks. If consecutive cases are considered, the next case is called for. 

COMMON BLOCKS 

COMMON blocks COFIIF, CSEVAL, INPUT, INTER, LOOKUP, 
OUTPUT, SAVED, and TABLES are used. 

TBL SUBROUTINES 

Subroutines BARCON, BARPRO, BARSET, CFEVAL, READIN, 
SEVAL, STARTS, and XNTERP call QUITS. 

QUITS calls subroutine DIRECT. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines or built-in FORTRAN functions are 
used. 

CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL QUITS 
SOLUTION METHOD _ 

Print out error statement. 

1. WRITE error message 

Print out COMMON block /COFJIF/. 

2. WRITE IFINT, AFINT, BFINT, CFINT, MMINT, TFINT 
Print out COMMON block /INPUT/. 
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;i. WRITE IDXMAX, ICTAB, IPRINT, ITWTAB, 1XTAB, MZETA, 

DXMAX, EPSZ, FJ, G, GAMO, PO, PHD, PIE, PRANDT, 
RBAR, SCALE, TO, THETA I, TOLCFA, TOLZET, 
TOLZME, ZMIIO, ZMVIS, ZNSTAN 

Print out COMMON block /OUTPUT/. 

4. WRITE B DELTA, CF, CI1, DELTA, DELSOT, DELSTR, FLAT, 
FORCE, IIG, PE, PHI, QW, SUMQDA, TE, THETA, 

TW, BE, X, XLARC, YR, Zl, Z2, Z3, Z4, Z5, ZETA, 
ZME 


Print out COMMON block /CSEVAL/. 

5. WRITE NOCTAB, IS, ROJ, FJG, CJG, GM102, GOGM1, POMAX, 

CPO, HO, SO, TCTAB, CP TAB, BCP, CCP, DCP, 

GTAB, II TAB, BARB1, BARB2, BARB 3 

Print out COMMON block /INTER/. 

6. WRITE CFAGT, CFAGP, CHPAR1, DX, DXRHO, HE, HW, 

IBEG, MZETAM, OOMZET, PHIP, PRE103, RHOE, 
RIIOUE, RMZETA, THETAP, XIBASE, XIEND, ZETATM, 
ZMZETA, ZMZETM, ZMZETP 

Print out COMMON block /SAVED/. 

7. WRITE A, B, C, ZI1, ZI1P, ZI2, ZI2P, ZI3, ZI3P, ZI4, ZI5, 

ZI6, ZI7 

Print out COMMON block /LOOKUP/. 

8. WRITE IMX, IPX, ITX, ITWX, IXPOS, IYX, CMX, CPX, CTX, 

CTWX, CYX 


Check whether Mach number tables have been input. 

9. If IX TAB 0, go to 34 

If IX TAB > 0, go to 10 

Check whether dimension of Mach tables supersedes the maximum 
value. 
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10. If IXTAB < 500, go to 15 
If IXTAB 2: 500, go to 11 

Sol number of table values to be printed, 

11. IS “ 495 

12. Go to 14 

Set number of table values to be printed to input table dimension. 

13. 13 - IXTAB 

Increase the number of table values to be printed. 

14. 13 = ( 10) ( (13) /( 10) + 1) 

Print out heading for tables printout. 

15. WRITE message for tables output 

Print out all the values in the axial distance table. 

16. Do 18, I - 1, 13, 10 
Compute printout counter. 

17. K = I + 9 

18. WRITE I, (XITAB(J), J = I, K) 

Print out all the values in the nozzle radius table. 

19. Do 21, I = 1, 13, 10 
Compute printout counter. 

20. K = I + 9 

21. WRITE I, (YITAB(J), J = I, K) 
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Print out all the values in the Mach number table. 

22. Do 24, I - 1, 13, 10 
Compute printout counter. 

23. K = I + 9 

24. WRITE I, (ZMTAB(J), J = I, K) 

Print out all the values in the wall temperature table. 

25. Do 27, I - 1, 13, 10 
Compute printout counter. 

26. K = I + 9 

27. WRITE I,. (TWTAB(J), J = I, K) 

Print out all the values of the calculated pressure table. 

28. Do 30, I = 1, 13, 10 

29. K = I + 9 

30. WRITE I, (PITAB(J), J = I, K) 

Print out all the values of the calculated temperature table. 

31. Do 33, I = 1, 13, 10 
Compute printout counter. 

32. K = I + 9 

33. WRITE I, (TITAB(J), J = I, K) 

Execute the next case, 

34. CALL DIRECT 

Subroutine QUITS nomenclature is given in Table 12. 
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TABLE 12. SUBROUTINE QUITS NOMENCLATURE 


Symbol 

A 

AFINT~ 


BA RBI 
"BARB 2 

BARB 3 


BDELTA 

BFINT 


Description 

Units 

Reference 

Wall enthalpy II 

w 

ft 2 /see 2 

/SAVED/, 7 

Wall enthalpy for 
integral evaluation 

ft 2 /sec 2 

/COFIIF/, 2 

Difference between 
stagnation enthalpy 
and wall enthalpy 

ft? /see 2 

/SAVED/, 7 

C polynomial 
equation 

Btu/ ( lbm-° R) 

/CSEVAL/, 5 

Negative value of 
the first derivative 
of BA RBI 

Etu/(lbm-° R 2 ) 

/CSEVAL/, 5 

Minus one -fourth the 
value of the second 
derivative of BA RBI 

Btu/(lbm-°R 3 ) 

/CSEVAL/, 5 

Coefficients in the 

Btu/ ( lbm-° R) 

/CSEVAL/, 5 


C - T relationship 

Boundary layer 
temperature thickness 
Difference between 
stagnation and wall 

enthalpy 

Difference between 
static and stagnation 
enthalpy 


ft 

ftVsec 2 


CCP 

Coefficient in the 
Cp - T relationship 

CF 

Skin friction 
coefficient 

CFAGP 

Skin friction 
coefficient obtained 
by using the Reynolds 
number R 

0 

CFAGT 

Initial value of skin 
friction coefficient 

CFINT 

Difference between 
static and stagnation 
enthalpy 



Btu/(lbm-®R 3 ) 



ft 2 /sec 2 


/OUTPUT/, 4 
/COFIIF/, 2 

/SAVED/, 7 


/CSEVAL/, 5 


/OUTPUT/, 4 


INTER/, 6 


/INTER/, 6 


/COFIIF/, 2 















TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Continued) 


Symbol 



CTWX 



Description 


Stanton number 


Term in the 
denominator of the 
Stanton number 

equation 

Specific heat at 
stagnation condition 

in work units 

Array of parabola 
coefficients for the 
Mach number table 

(ZMTAB) 

Specific heat 
coefficient at constant 
pressure for stagna- 
tion condition 

Array of C values 
P 

corresponding to the 
values in the temper- 
ature table 

Array of parabola 
coefficients for the 
contour point pres- 
sure table (PITAB) 
Array of parabola 
coefficients for the 
wall temperature 
table (TWTAB) 

Array of parabola 
coefficients for the 
contour point temper- 
ature table (TITAB) 
Array of parabola 
coefficients for the 
nozzle radius table 
(YITAB) 


Units 



ft 2 /(sec 2 - 0 R) 


lbm-°R) 


Btu/ ( lbm-° R) 


Reference 


/OUTPUT/, 4 


/INTER/,' 6 


/CSEVAL/, 5 


/LOOKUP/, 8 


/CSEVAL/, 5 


/CSEVAL/, 5 


/LOOKUP/, 8 


/LOOKUP/, 8 


/LOOKUP/, 8 


'LOOKUP/, 8 





















TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Continued) 


Symbol 


DCP 


DELSOT 


DELSTR 


DELTA 


DX 


DXMAX 


DXRHO 


EPSZ 

FJ 


FJG 


FLAT 


FORCE 


GAMO 


GM102 


Description 


Coefficient in the 
C - V relationship 


Boundary layer 
shape factor 


Boundary layer 

displacement 

thickness 


Boundary layer 
velocity thickness 


Weighted difference 
of tabic values of the 
axial distance 


Maximum length of 
step size 


One-tenth the 
difference between 
axial distance points 


Flow geometry 
indicator 


Conversion factor 
between thermal 
and work units 


FJ multiplied by G 


Force normal to x- 
direction for two- 
dimensional planar 
flow 


Drag force in axial 
or x-direction 


Acceleration of 
gravity used as a 
proportionality 
constant 


Specific heat ratio at 
stagnation condition 


One-half the specific 
heat ratio at stag- 
nation condition 
minus one 


Units 


Btu/( Ibm J* R) 


ft 


ft 


ft 


ft 


ft 


( ft-lbf) /Btu 


( ft 2 -lbm) / ( Btu-s ec 2 ) 


lbf 


lbf 


( lbm-ft) / ( lbf-sec 2 ) 


Reference 


/CSEVAL/, 


/OUTPUT/, 4 
"/OUTPUT/, 4 


/OUTPUT/, 4 


/INTER/, G 


/INPUT/, 3 


/INTER/, 6 


/INPUT/, 3 


/INPUT/, 3 


/CSEVAL/, 


/OUTPUT/, 4 


/OUTPUT/, 4 


/INPUT/, 3 


/INPUT/, 3 


/CSEVAL/, 5 


Vw 
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TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Continued) 


Symbol 

Description 

Units 


GOG Ml 

Specific heat ratio at 
stagnation condition 
divided by the specific 
heat ratio minus one 

— 


GTAD 

Array of terms used 
by SEVAL. 

Btu/ ( lbm-° It) 

/CSEVAL/, 5 

no 

Stagnation 

enthalpy 

ft 2 /sec 2 

/CSEVAL/, 5 

1IE 

Enthalpy at static 
temperature 

ftVsec 2 

/INTER/, g 

IIG 

Heat transfer 
coefficient 

Btu/( ft 2 -sec-° R) 

/OUTPUT/, 4 

IITAB 

Array of enthalpy 
values 

Btu/lbm 

/CSEVAL/, 5 

HW 

Wall enthalpy at wall 
temperature 

ft 2 /sec 2 

/INTER/, 6 

I 

Do-loop counter 

— 

16-33 

13 

Number of table 
values to be printed 
out 


11, 13, 14, 
16, 19, 

22, 25, 28, 
31 

IBEG 1 

Subscript counter at 
which Mach number 
table exceeds one 


/INTER/, 6 

ICTAB 

Specific heat indica- 
tor or table dimension 

— 

/INPUT/, 3 

IDXMAX 

Step-size indicator 

— 

/INPUT/, 3 

IFINT 

Integral evaluation 
indicator 

— 1 

/COFIIF/, 2 

I MX 

Position or start 
indicator for Mach 
number table 


/LOOKUP/, 8 

IPRINT 

Printout indicator 

— 

/INPUT/, 3 

IPX 

Position or start 
indicator for contour 
point pressure table 


/LOOKUP/, 8 

IS 

One less than the 

dimension of the 

C - T 
P 


/CSEVAL/, 5 
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TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

ITWTAB 

Wall temperature 
indicator 

— 

/INPUT/, 3 

ITWX 

Position or start 
indicator for wall 
temperature tabic 

1 

/LOOKUP/, 8 

ITX 

Position or start 
indicator for contour 
point temperature 
table 


/LOOKUP/, 8 

IXPOS 

Array position 
indicator 

— 

/LOOKUP/, 8 

KTAB 

Dimension of x, y, 
Mach tables 

— 

/INPUT/, 3 

IYX 

Position or start 
indicator for nozzle 
radius table 

~ 

/LOOKUP/, 8 

J 

Subscript printout 
counter 

— 

18, 21, 24, 
27, 30, 33 

K 

Printout counter 


17, 18, 20, 
21, 23, 24, 
26, 27, 29, 
30, 32, 33 

MMINT 

Exponent for integral 
evaluation 

— 

/COFIIF/, 2 

MZETA 

Exponent of 
velocity profile 

— 

/INPUT/, 3 

MZETAM 

MZETA minus one 

— 

~ /INTER/, 6 J 

NOCTAB 

Saved value of 
ICTAB 

— 

" /CSEVAL/, 5 

OOMZET 

One divided by 
Z MZETA 

— 

/INTER/, 6 

PO 

Stagnation 

pressure 

lbf/ft 2 

/INPUT/, 3 

POMAX 

Saved value of 
stagnation pressure 

lbf/ft 2 

/CSEVAL/, 5 

PE 

Static pressure 

Ibf/ft 2 

~ /OUTPUT/, 
4 

PHI 

Boundary layer 
energy thickness 

ft 

/OUTPUT/, 4 

PHII 

Input or initial value 
of energy thickness 

ft 

7INPUT/, 3 
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TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Continued) 


Symbol 


PHIP 


PIE 


PITAB 


PRANDT 
PRE103 ' 


QW 


RBAR 


Description 


Slope of energy 
thickness 


Circular constant tt 


Units 


Array of pressures 
at contour points 
x. y 


Prandtl number 


Recovery factor 


Local heat transfer 
rate to the wall 


Specific gas 
constant 


lbf/ft 5 


Btu/ ( ft 2 -sec) 


( ft-lbf) / ( lbm-° R) 


Reference 


/INTER/, G 
/INPUT/. 3~ 


/TABLES/, 30 


/INPUT/, 3 
/INTER/. 6 


/OUTPUT/, 4 


/INPUT/, 3 


RHOE 


RHOUE 


Density at Boundary 
layer c-dge 


lbm/ft 3 


Flow density (p • v) 
at boundary layer 
edge 


lbm/(ft z -scc) 


/INTER/, 6 


/INTER/, 6 


RMZETA 


One divided by 
ZMZETP 


/INTER/, 6 
/CSEVAL/ , 5 


ROJ 


SO 


Specific gas 
constant 


Btu/(lbm-°R) 


Stagnation entropy 


Btu/ ( lbm-° R) 


/CSEVAL/, 5 


SCALE 


Contour scale factor 


or ft 


/INPUT/, 3 


'/OUTPUT/, 4 


SUMQDA 


Integrated heat 
transfer rate 


Btu/sec 


TO 


TCTAB 


Stagnation 

temperature 


R 


Temperature table 

for C values 

P 


’R 


/INPUT/, 3 
/CSEVAL/7T 


TTT 


TFINT 


Static temperature 


’R 


Static temperature 
for integral evaluation 


R 


/OUTPUT/~ 


/COFIIF/, 2 


THETA 


THETA 


Boundary layer 
momentum thickness 


ft 


Boundary layer 
momentum thickness 


ft 


/OUTPUT/, 4 


/OUTPUT/, 4 
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TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Continued) 


Symbol 


THETA I 


TIIETAP 


TITAB 


TOLCFA 


TOLZET 


TOLZME 


TW 


TWTAB 


UE 


X 


XIBASE 


XIEND 


XITAB 


XLARC 


YITAB 


YITAB 


Description 

Input or initial 
value of momentum 

thickness 

Slope of momentum 

thickness 

Array of tempera- 
tures at contour 

points, x t y 

Tolerance in 

C f - C f R s 
iteration 


Tolerance in the 
zeta interation 


Units 



Wall temperature 


Array of wall tem- 
peratures at contour 
points x, y 


Velocity of fluid 


Axial distance 


First value of 
axial distance 
table 


Last value of axial 
distance table 


Array of axial 
distance values 


Arc length of the 

contour 

Array of nozzle 
radius or contour 

height 

Array of nozzle 
radius or contour 
height 


Reference 


/INPUT/, 8 


/INTER/, 6 

/TABLES/, 

33 

/INPUT/, 3 



/OUTPUT/, 4 


/TABLES/, 

27 



TABLES/, 21 


/TABLES/, 

21 














































TABLE 12. SUBROUTINE QUITS NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

ZI6 

I. - j* — S n els 
o 

— 

/SAVED/, 7 

ZI7 

s n ds 

£ p o 

— 

/SAVED/, 7 

ZME 

Mach number 



/OUTPUT/, 4 

ZMTAB 

Array of Mach numbers at the 
contour points, x, y 

— 

~ /TABLES/, 
24 

ZMUO 

Viscosity at stagnation 
temperature 

mg! 


ZMVIS 

Exponent in viscosity- 
temperature law 


/INPUT/, 3 

ZMZETA 

Real value of MZETA 

— 

/INTER/, 6 

ZMZETM 

ZMZETA minus one 

. .. 

/INTER/, 6 

ZMZETP 

ZMZETA plus one 

— 

/INTER/, 6 

ZNSTAN 

Interaction exponent in 
Stanton number 
relationship 


/INPUT/, 3 














SUBROUTINE READIN 


Subroutine READIN first calls for single item variables such as 
integers, constants, and gas constants. It then reads in tabular values such 
as contour points, Mach number, and wall temperature. It prints out the 
input and performs diagnostics to detect obvious inconsistencies in input 
quantities. 

Before any input is read in, the following constants are set to their 
nominal values and need be input only if different values are desired. 


1. 

Contour scale factor 

SCALE =1.0 

2. 

Velocity power law 
exponent 

MZETA - 7 

3. 

Interaction exponent in 
Stanton number relation 

ZNSTAN = 0. 10 

4. 

Conversion factor between 
thermal and work units 

FJ = 778.20 (ft-lb)/Btu 

5. 

Sea level acceleration of 
gravity 

G = 32.1740 (ft-lbm)/(lbf-sec 2 ) 

6. 

Tolerance in skin friction 
coefficient iteration loop 

TOLCFA ■ 1.0 x 10" 4 

7. 

Tolerance in gas property 
evaluation loops 

TOLZME ■ 1.0 X 10 -7 

8. 

Tolerance in zeta iteration 
loop 

TOLZET = 0.00030 

The first input card for each case is a title card containing Hollerith 
information in columns 1 through 78. Input data are read in by the NAMELIST 
name NAM1. The variables that are included in the NAM1 NAMELIST are 
as follows: CPTAB, DXMAX, EPSZ, FJ, G, GAMO, ICTAB, IPRINT, 
ITWTAB, IXTAB, MZETA, PO, PHII, PRANDT, RBAR, SCALE, TO, 
TCTAB, THETAI, TOLCFA, TOLZET, TOLZME, TWTAB, XITAB, 
ZMTAB, ZMUO, ZMVIS, ZNSTAN 


COMMON BLOCKS 




COMMON blocks CSEVAL, INPUT, and TABLES are used. 

TBL SUBROUTINES 

Subroutine DIRECT calls READIN. 

REA DIN calls subroutine QUITS. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines are used. 

Built-in FORTRAN function IABS is used, 

CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL READIN 

SOLUTION METHOD 

Set constants to nominal values. 

1. SCALE - 1.0 

2. MZETA = 7 

3. ZNSTAN = 0. 10 

4. FJ = 778.20 

5. G = 32.1740 

6. TOLCFA = 1.0 X 10“ 4 

7. TOLZME = 1.0 X 10" 7 

8. TOLZET = 0.00030 

Save value of maximum step size from previous case. 
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9. DXMAXO - DXMAX 


Set maximum stop size to check whether saved value from 
previous case is to be used. 

10. DXMAX = -28982.0 

Head in title for this case. 


11. HEAD TITLE 

Read in name list data for this case, 

12. READ ($NAMl) 

Check whether a new value for the maximum step size has been input. 

13. If DXMAX + -28982.0, go to 17 
If DXMAX = -28982.0, go to 14 

Check whether step site from previous case or computed value must be 
used. 

14. If IDXMAX = 0, go to 21 
If IDXMAX * 0, go to 15 

Set maximum step size to saved value from previous case. 

15. DXMAX = DXMAXO 

16. Go to 22 

Step size has been input; cheek whether actual step size must be 
computed. 


17. If DXMAX ^ 0.0, go to 20 
If DXMAX > 0.0, go to 18 
Set step size indicator to use the same step size. 
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18. IDXMAX - 1 


19. Go to 22 

Set step size indicator to use a new step size. 

20. IDXMAX = 0 

Compute value of maximum step size to be used. 

21. DXMAX = ((XITAB(IXTAB) - XITAB( l) )/( 100. )) (SCALE) 
Write out the title heading for this case. 

22. WRITE TITLE 

Initialize error indicator to zero. 

23. IERROR = 0 

Print out the velocity profile exponent. 

24. WRITE MZETA 

Check value of velocity profile exponent. 

25. If MZETA 2= 0, go to 28 
If MZETA < 0, go to 26 

If the velocity profile exponent is in error, print out message. 

26. WRITE error message about input value 
Set input error indicator to one. 

27. IERROR = 1 

Write out printout indicator. 

28. WRITE IPRINT 
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Chock value of printout indicator. 




29. If IP HINT = 1 or IPUINT = 0, go to M2 
If IPUINT t- 1 and IPUINT * 0, go to 20 
The printout indicator is in error; print out message. 

90. WRITE error message about input value 
Set input error indicator to one. 

31. IKUUOU = 1 

Print out dimension of x, y, and M tables. 

32. WRITE IXTAB 

Check whether value of dimension indicator is within limits. 

33. If IXTAB 2= 4 and IXTAB ^ 500, go to 36 
If IXTAB < 4 or IXTAB > 500, go to 34 

Dimension indicator is below minimum limit or above maximum limit. 

34. WRITE error message about input value 
Set input error indicator to one. 

35. IERROR = 1 

Print out specific heat indicator. 

36. WRITE ICTAB 

Check specific heat option to be used. 

37. If ICTAB = 0, go to 41 
If ICTAB * 0, go to 38 



Tabular specific heat option is called for, check whether dimension 
is within limits. 

38. If ICTAB 3 and ICTAB r= 20, go to 41 
If ICTAB < 3 or ICTAB > 20, go to 30 

Dimension indicator is out of limits. 

39. WRITE error message about input value 
Set input error indicator to one. 

40. IERROR = 1 

Write out wall temperature option indicator. 

41. WRITE ITWTAB 

Check on value of wall temperature option indicator. 

42. If | ITWTAB | = 1 or ITWTAB = 0, go to 45 
If I ITWTAB | * 1 and ITWTAB * 0, go to 43 

Wall temperature indicator is in error. 

43. WRITE error message about input value 
Set input error indicator to one. 

44. IERROR = 1 

Print out the value of the stagnation temperature. 

45. WRITE TO 

Check value of stagnation temperature for error. 

46. If TO > 0.0, go to 49 
If TO s o.O, go to 47 




Stagnation temperature value is in error, 

47. WRITE error message about input value 
Set input error indicator to one, 

48. IKRROR = 1 

Print out the value of the stagnation pressure. 

49. WRITE PO 

Check value of stagnation pressure for error. 

50. If PO > 0.0, go to 53 
If PO =£ 0.0, go to 51 

Stagnation pressure value is in error. 

51. WRITE error message about input value 
Set input error indicator to one. 

52. IERROR = 1 

Print out stagnation value of the specific heat ratio. 

53. WRITE GAM0 

Check on specific heat ratio option indicator, 

54. If ICTAB * 0, go to 58 
If ICTAB = 0, go to 55 

Constant specific heat is to be used, check value of stagnation specific 
heat ratio. 

55. If GAM0 > 1.0, go to 58 


If GAM0 l.o, go to 56 


Stagnation specific heat ratio is in error. 

56 . WRITK error message about input value 
Set input error indicator to one. 

57 . IKRHOR = 1 

Print out the value of the Prandtl number. 

58. WRITE PRANDT 

Check on value of the Prandtl number. 

59 . If PRANDT > 0.0, go to 02 
If PRANDT * 0.0, go to 00 

The Prandtl number is in error. 

60. WRITE error message about input value 

Set input error indicator to one. 

61. IERROR * 1 

Print out the stagnation viscosity. 

62. WRITE ZN 00 

Check on value of the stagnation viscosity. 

63. If ZMU0 > 0.0, go to 66 
If ZMU0 ^ 0.0, go to 64 

The stagnation viscosity is in error. 

64. WRITE error message about input value 


Set input error indicator to one. 


65. 1ERR0R = 1 


Print out the temperature-viscosity relation exponent. 

66. WRITE ZMVIS 

Print out the Stanton number interaction exponent. 

67. WRITE ZNSTAN 

Print out the maximum length of the step size. 

68. WRITE DXMAX 

Check whether the sonic point start procedure is to be used. 

69. If THETAI < 0. 0, go to 72 . 

If THETAI ^ 0.0, go to 70 

Print out the initial value of the momentum thickness, 

70. WRITE THETAI 

Print out the initial value of the energy thickness. 

71. WRITE PHII 

Print out the flow geometry indicator. 

72. WRITE EPSZ 

Check value of the flow geometry indicator. 

73. If EPSZ = 0.0 or EPSZ - 1.0, go to 76 
If EPSZ * 0.0 and EPZS * 1.0, goto 74 

The flow geometry indicator is in error. 

74. WRITE error message about input value 
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Set input error indicator to one. 


75. IERROR = 1 

Print out the value of the specific gas constant. 

76. WRITE REAR 

Check value of specific gas constant. 

77. If RBAR > 0.0, go to 80 
If RBAR =s 0.0, go to 78 

The specific gas constant is in error. 

78. WRITE error message about input value 
Set input error indicator to one. 

79. IERROR = 1 

Print out the conversion factor between thermal and work units. 

80. WRITE FJ 

Check value of the conversion factor. 

81. If FJ > 0.0, go to 84 
If FJ £ 0. 0, go to 82 

The conversion factor is in error. 

82. WRITE error message about input value 
Set input error indicator to one. 

83. IERROR = 1 

Print out the value of the acceleration of gravity at sea level. 

84. WRITE G 
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Check value of the acceleration of gravity. 

85. If G > 0.0, go to 88 
If G s 0. 0, go to 86 

The value of the acceleration of gravity is in error. 

86. WRITE error message about input value 
Set input error indicator to one, 

87. IERROR = 1 

Print out the nozzle contour scale factor. 

88. WRITE SCALE 

Check whether the skin friction coefficient tolerance has been input. 

89. If TOLCFA = 1.0 x 10 -4 , go to 91 
If TOLCFA # 1.0 x 10~ 4 , go to 90 

Print 'out the input value of the skin-f.niction coefficient tolerance. 

90. WRITE TOLCFA 

Check whether the zeta iteration tolerance has been input. 

91. If TOLZET = 0.00030, go to 93 
If TOLZET * 0.00030, goto 92 

Print out the input value of the zeta iteration tolerance. 

92. WRITE TOLZET 

Check whether the gas property evaluation tolerance has been input. 

93. IF TOLZME * 1.0 x 10 -T , go to 95 
IF TOLZME # i.o x 10~ T , go to 94 
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Print out the input value of the gas property evaluation tolerance. 

94. WRITE TOIiZME 

Check whether the specific heat table has been input. 

95. If ICTAB =s o, go to 112 
If ICTAB > 0, gotoM 

Print out a heading for the specific heat/temperature tables. 

96. WRITE heading for specific heat/temperature tables 
Print out the specific heat/temperature tables. 

97. WRITE (I, CPTAB(l) , TCTAB(l), 1= 1, ICTAB) 

Set limit for do-loop. 

98. n = ICTAB - 1 

99. Do 103, I = 1, II 

Check whether temperature table values are monatonically increasing. 

100. If TCTAB (I + 1) > TCTAB (I), goto 103 
If TCTAB (I + 1) =£ TCTAB (I), goto 101 

Temperature table values are not in the proper order. 

101. WRITE error message about temperature table 
Set input error indicator to one. 

102. IERROR = 1 

103. Continue 

Check the first value in the temperature table. 
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104. If TCTAB (1) > 0.0, go to 107 
If TCTAB (1) S 0.0, go to 105 

The first temperature table value is in error. 

105. WRITE error message about temperature table 
Set input error indicator to one, 

106. I ERROR = 1 

107. Do 111, I = 1, ICTAB 

Check whether specific heat table values are greater than zero. 

108. If CPTAB (I) > 0.0, goto 111 
If CPTAB (I) =£ 0.0, go to 109 

Specific heat table values are in error. 

109. WRITE error message about specific heat table 
Set input error indicator to one. 

110. IERROR = 1 

111. Continue 

112. Do 116, I = 1, KTAB 
Check Mach number table values. 

113. If ZMTAB (I) > 0.0, go to 116 
If ZMTAB (I) ^ 0.0, go to 114 

The Mach number table is in error. 

114. WRITE error message about Mach number table 
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Set input error indicator to one. 


115. I ERROR = 1 

116. Continue 

Set limit for do-loop. 

117. II = IXTAB - 1 

118. Do 122, I = 1, II 

Check whether axial distance table values are monatonically increasing. 

119. If XITAB (I + 1) & XITAB (i), goto 122 
If XITAB (I + 1) < XITAB (i), goto 120 

The axial distance table values are not in proper order. 

120. WRITE error message about the axial distance table 

Set input error indicator to one. _ 

121. IERROR = 1 

122. Continue 

Check wall temperature option indicator. 

123. If ITWTAB < 0, go to 133 
If ITWTAB = 0, go to 130 
If ITWTAB > 0, go to 124 

Tabular wall temperatures are to be used. 

124. Do 128, I = 1, IXTAB 

Check whether wall temperature values are greater than zero. 
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125. If TWTAB (I) > 0.0, goto 128 
If TWTAB (I) £ 0.0, goto 126 

The wall temperature table is in error. 

126. WRITE error message about wall temperature table 
Set input error indicator to one, 

127. IERROR = 1 

128. Continue 

129. Go to 133 

Constant wall temperature option is to be us£d; check for proper 
values. 

130. IfTWTAR(l) > 0.0, go to 133 
IfTWTAB(l) < 0.0, go to 131 

Input value of constant wall temperature is in error. 

131. WRITE error message about wall temperature 
Set input error indicator to one. 

132. IERROR = 1 

Check whether contour scale factor has been input. 

133. If SCALE = 1.0, go to 137 
If SCALE * 1.0, go to 134 

Scale factor has been input} calculate real values for x and y. 

134. Do 136, I = 1, IXTAB 

135. XITAB(I) = (XITAB(I)) (SCALE) 
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136. YITAB (I) = (YITAB(I))(SCALE) 

Check wall temperature option indicator. 

137. If ITWTAR < 0, go to 139 
If ITWTAR = 0, go to 138 
If ITWTAB > 0, go to 142 

Print out value of constant wall temperature to be used. 

138. WRITE TWTAB(l) 

Print out heading for x, y, and M tables. 

139. WRITE heading for x, y, and M tables 
Print out input values of x, y, and M tables. 

140. WRITE (I, XITAB(I), YITAB(I), ZMTAB(I) , I - 1. IXTAB) 

141. Go to 144 

Print out heading for x, y, M, and T^ tables. 

142. WRITE heading for x, y, M, and T tables 

w 

Print out input values of x, y, M, and T tables. 

w 

143. WRITE (I, XITABtl), YITAB(I), ZMTAB(I) , 

TWTAB(I) ,1=1, IXTAB) 

Check dimension indicator for x, y, and M tables. 

144. If IXTAB 1, go to 150 
If IXTAB > 1, go to 145 

Check whether axial distance values are in monatonically increasing 
order. 

145. Do 149, I = 2, IXTAB 
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146. If XITAB (I) > XITAB (I - l) , go to 149 
If XITAB (I) ^ XITAB (I - l) , goto 147 

The axial distance table is in error. 

147. WRITE error message about axial distance table 
Set input error indicator to one. 

148. IERROR = 1 

149. Continue 

Check wall temperature option indicator. 

150. If ITWTAB < 0, go to 154 
If ITWTAB a 0, go to 151 

Check for sonic point start procedure option. 

151. If THETAI 2= 0.0, goto 154 
If THETAI < 0.0, goto 152 

The sonic point start option is incompatible with constant or tabular 
wall temperature options. 

152. WRITE error message about input options 
Set input error indicator to one. 

153. IERROR = 1 

Check for any errors in the input values. 

154. If IERROR ^ o, return 
If IERROR > 0, go to 155 

There is an error in the input; print out COMMON blocks, and go to 
the next case. 

155. CALL QUITS 
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Subroutine READIN nomenclature is given in Table 13. 


TABLE 13. SUBROUTINE READIN NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

CPTAB 

Array of specific 
heat values 

Btu/( lbm-° R) 

/CSEVAL/, 
/NAM1./, 
97, 108 

DXMAX 

Maximum length of 
step size 

ft 

/INPUT/, 
/NAMI /, 

9, 10, 13, 
15, 17, 21, 
68 

DXMAXO 

Maximum step size 
from last case 

ft 

9, 15 

EPSZ 

Flow geometry 
indicator 


/INPUT/, 
/NAMI/, 
72, 73 

FJ 

Conversion factor 
between thermal 
and work units 

(ft-lbf) /Btu 

/INPUT/, 
/NAMI/, 
4, 80, 81 

G 

Acceleration of 
gravity at sea level 

( ft-lbm) /( lbf-sec 2 ) 

/INPUT/, 
/NAMI/, 
5, 84, 85 

GAMO 

Specific heat ratio 
at stagnation 
condition 


/INPUT/, ' 

/NAMI/, 

53, 55 

I 

Do-loop and print- 
out counter 


99, 100, 107, 
108, 112, 

113, 118, 

119, 124, 

125, 134- 
136, 140, 

143, 145, 

146 

ICTAB 

Specific heat 
indicator and/or 
table dimension 


/INPUT/, 
/NAMI/, 
36-38, 54, 
95, 97, 98, 
107 
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TABLE 13. SUBROUTINE READIN NOMENCLATURE ( Continued) 


Symbol 

Description 

Units 

Holrrorn’o 

IDXMAX 

Step size indicator 

— 

/INPUT/, 

14, 18, 20 

23, 27, 81 , 

35, '10, '14, 

48, !>2, 57, 

(11, 05, 71), 

70, 88, 87, 

102, 100, 

11(1, 115, 

121, 127, 

182, 148, 

153, 154 

I ERROR 

Input error 
indicator 


II 

Do-loop limit 

— 

98, 99, 
117, 118 

IPRINT 

Print -out 
indicator 


/INPUT?, 
/NAM1/., 
28, 29 

ITWTAB 

Wall temperature 
option indicator 1 


/INPUT/, 

/NAM1/, 

41, 42, 123, 
137, 150 

IXTAB 

Dimension indicator 
for x, y, and M 
tables 


/INPUT/, 

/NAM1/, 

32, 33, 112, 
117, 124, 
134, 140, 
143-145 

MZETA 

Velocity profile 
exponent 


/INPUT/, 
/NAMl/, 2, 
24, 25 

PO 

Stagnation pressure 

lbf/ft 2 

/INPUT/, 
/NAMl/, 
49, 50 

PHII 

Initial .value of 
energy thickness 

ft 

/INPUT/, 
/NAMl/, 71 

PRANDT 

Prandtl number 


/input/, 

/NAMl/, 
58, 59 
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Symbol 


Description 


Units 


Reference 


RBAR 


SCALE 


TO 


TCTAB 


Specific gas 
constant 


Contour scale 
factor 


Stagnation 

temperature 

Array of tem- 
perature values 
corresponding to 
the specific heat 


( ft-lbf) / ( lbm-° K) 


°R 


°R 


/INPUT/, 

/NAM1/, 

76. 77 
/INPUT/, 
/NAM1/, 1, 
21, 88, 133, 
135, 136 
/INPUT/, 
/NAM1/, 

45, 46 

/CSEVAL/, 

/NAM1/, 

97, 100, 104 


THETAI 


table values 

Initial value of 
momentum thickness 


ft 


TITLE 


Array of Hollerith 


/INPUT/, 

/NAM1/, 

69, 70, 151 
DIM, 11, 22 


characters input 
as the title heading 
for the case being 


TOLCFA 


run 

Skin friction 
iteration tolerance 


TOLZET 


/INPUT/, 
/NAM1/, 
6, 89, 90 


Zeta iteration 
tolerance 


TOLZME | Gas property 

evaluation tolerance 


/input/, 

/NAM1/, 8, 
91, 92 


/input/, 

/NAM1/, 7, 
93, 94 


TWTAB | Array of wall 

temperature values 


/TABLES/, 
/NAMi/, 
125, 130, 
138, 143 


t 


i 
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TABLK 13. SUBROUTINE READIN NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

XITAB 

Array of axial 
distance values 

ft 

/TABLES/, 
/NAM1/, 
119, 135, 
140, 143, 
146 

YITAB 

Array of nozzle 
radius values 

ft 

/TABLES/, 

/NAMl/, 

136, 140, ; 

143 

ZMTAB 

Array of Mach 
number values 


/TABLES/, 
/NAMl/, 
113, 140, 
143 

ZMUO 

Stagnation 

viscosity 

lb/(ft,-sec) 

/INPUT/, 
/NAMl/, 
62, 63 

ZMVIS 

Viscosity- 
temperature law 
exponent 

1 

/INPUT/, 

/NAMl/, 

66 

ZNSTAN 

Stanton number 
interaction exponent 


/INPUT/, 

/NAMl/, 

67 


SUBROUTINE SEVA I. 


Subroutine SEVA 1 j provides the evaluation of several thermodynamic 
properties under the basic assumptions that the enthalpy is a function of 
temperature (thermally perfect) only and the C is expressed as an analytic 

function between intervals of a C - T table, the function being a piecewise 
cubic. 15 


The basic relations connecting the thermodynamic quantities are: 

C p = C pi + DCP i ( T ~ V + CCP i (T " T i^ 2 + DCP i ( T - T P 3 * (!) 


T 

H (T) = H + / C (t)dt , (2) 

T 

i 

T C (t) 

S(T, P) = S(T P ) + / -£■ — dt - H« In ±- . ( 3 ) 

° o T C o 

o 

These functions are solved for several modes of input to SEVAL by an 
option indicator IND1. 


JND1 - 0 S is evaluated directly from equation (3) for given values 
of P and T . 

S= S(T, P) . 

• SEVAL (0, T, P, S) 

IND1 = 1 H is evaluated directly from equation (2) for a given 
T, and is evaluated directly from equation (1) . 

H = H( T) , C = C ( T) . 

P P 

• SEVAL (1. T. CP. II) 


IND1 = 2 T is evaluated from equation (2) for a given H in an 
inverse manner using a modified pseudo-position 
procedure. Then C is evaluated directly from 
equation ( 1) . p 
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T = T(H), C = C (T) . 

P P 

• SEVAL (2, T, CP, H) 

IND1 =3 T is evaluated from equation (3), for given S and P, in 
an inverse manner using a modified pseudo-position 
procedure. 


T - T (S, P) 

• SEVAL (3, T, P, S) 

IND1 = -1 P is evaluated from equation (3) written in the form 


In 




C (t) 

-E-r— dt + S(T , P ) - S(T, P) 
t oo 


/R* 


P = P (T, s) 

• SEVAL (-1, T, P, S) 


COMMON BLOCKS 

COMMON blocks CSEVAL and INPUT are used. 

TBL SUBROUTINES 

Subroutines BARPRO, BARSET, FIIF, GETPT, and START call 
SEVAL. 

SEVAL calls subroutine QUITS. 

FORTRAN SYSTEM ROUTINES 

FORTRAN library routines ALOG and EXP are used. 

Built-in FORTRAN function ABS is used. 

CALLING SEQUENCE 


The subroutine calling sequence is 
CALL SEVAL (IND1, AA, BB, CC) 


159 


where 


IND1 = option indicator (0, 1, 2, 3 or -l), 

AA = temperature T, 

BB = pressure P or specific heat C^, 

CC = entropy S or enthalpy H. 

Calling sequence depends on the option indicator represented by the 
first term in the parentheses, 

SEVAL (0, T, P, _S ) 

SEVAL (1, T, CP, H ) 

SEVAL (2, T, CP, H ) 

SEVAL (3, T, P, S ) 

SEVAL (-1, T, P, S ) 

SOLUTION METHOD 

Define the function of entropy at temperature T and pressure P. 

1. GAPF( T, G, Al, Bi, Cl, Dl) = G + (Al)(Ln(T)) + (B1)(T) 

+ (C1)(T) 2 + (Dl) (T) 3 / (3. 0) - PR 

Set three input values as 

2. T = AA 

3. A = BB 

4. B = CC 

Check the option indicator IND1. 
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5. If IND1 < 2, go to 11 
If IND1 = 2, go to 6 
If IND1 > 2 , go to 16 

Save the input enthalpy in thermal units 
(IND1 = 2; H -*■ T, CP-) 

6. B = (B)/(FJG) 

Check whether variable or constant specific heat is considered. 

7. If ICTAB > 0, go to 55 
If ICTAB =£ 0, go to 8 

Save temperature and specific heat in the case of constant specific 

heat ( C = C ) . 

P po 

8. T = (B)/(CP0) 

9. A = CPO 

10. Go to 106 

Check the option indicator IND1, 

11. If IND1 < 1, go to 33 
If IND1 s l, go to 12 

Check whether variable or constant specific heat (IND1 = 1; T C^> H) 
is used. 

12. If ICTAB > 0, go to 33 
If ICTAB =£ 0, go to 13 

Save specific heat. 

13. A = CPO i ^ 

1 


Calculate enthalpy in work units: (C gJT ) . 

po e 

14. B = (CJG)lT) 

15. Go to 106 

Calculate the last term in the entropy equation (thermal units). 

~ In (IND1 = 3; P, S - T) . 

16 . PR = ( ROJ) ( Ln( ( A)/ ( POMAX) ) ) 

Calculate entropy to search the temperature table. 

17. STAB = GAPF( TCTAB( IS) , GTAB(IS), BARBl(IS), 

BARB2US), BARB3( IS) , DCP(IS)) 

Check the values between input entropy and calculated entropy above. 

18. If B ^ STAB, go to 22 
If B < STAB, go to 19 

Decrement subscript counter to use next lower table values. 

19. IS = IS - 1 

Check whether subscript counter is still greater than zero. 

20. If IS ^ 0, go to 42 
If IS > 0, go to 21 

21. Go to 17 

Calculate the entropy corresponding to TCTAB(IS + l). 

22. STAB = GAPF( TCTAB( IS + 1), GTAB(IS) , BARBl(IS), 

BARB2US), BARB3( IS) , DCP(IS)) 


Check whether the calculated entropy above exceeds the input entropy 
value. 

23. If B < STAB, go to 27 
If B s STAB, go to 24 

Increment the subscript counter to use the next higher table values. 

24. IS = IS + 1 

Check whether the subscript counter has exceeded the table dimension. 

25. If IS s 20, go to 42 
If IS < 20, go to 26 

Repeat the above search until B < STAB is obtained. 

26. Go to 22 
Check the value of IS. 

27. If IS 5: NOCTAB or IS =£ 0, go to 42 
If IS < NOCTAB and IS > 0, go to 28 

Save values of temperature and entropy. 

28. TTP = TCTAB(IS) 

29. FP = GAPF(TCTABUS), GTAB(IS), BARBl(IS), 

BARB2(IS), BARB3( IS) , DCP(IS)) 

30. TTPP = TCTABdS + 1) 

31. FPP = STAB 

32. Go to 68 

a 
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M 


Temperature T is checked in the following steps down to 44, 

(IND1 = 0 or l). 

33. If T ^ TCTAB(IS), go to 37 
If T < TCTAB(IS) , go to 34 

34. IS = IS - 1 

35. If IS 0, go to 42 
If IS > 0, go to 36 

36. Go to 33 

37. If T < TCTAB(IS + l), go to 41 
If T ^ TCTAB(IS + 1), goto 38 

38. IS = IS + 1 

39. If IS & 20, go to 42 
If IS < 20, go to 40 

40. Go to 37 

41. If IS > 0, go to 44 
If IS ^ 0, go to 42 

Print errors and corresponding values, 

42. WRITE IS, IND1, T, A, B 

Stop calculation by calling subroutine QUITS. 

43. CALL QUITS 

Check whether IS exceeds-NOCTAB. 

44. If IS s NOCTAB, go to 42 
If IS < NOCTAB, go to 45 


Check the indicator for option IND1. 


45. If IND1 < 0, go to 53 
If IND1 = 0, go to 50 
If IND1 > 0, go to 46 

Set T minus tabulated value TCTAB.(IS) 

(IND1 = 1; T - CP, H) • 

46. DELT = T - TCTAB(IS) 

Calculate enthalpy in thermal units according to the equation in 
subroutine BARSET. 

47. B = HTABt IS) + (CPTAB( IS) ) ( DELT) 

+ ( 0. 5) ( BCP( IS)) ( DELT) 2 + (CCP(IS)) (DELT) 3 /(3. 0) 

+ (0. 25) (DCP(IS)) (DELT) 4 

Save enthalpy in work units. 

48. B = ( B) ( FJG) 

49. Go to 105 

Calculate entropy in this case (IND1 = 0; T, P — S) . 

50. PR = ( ROJ) ( Ln( l A)/ P0MAX) ) ) 

Entr opy: 

51. B = GAPF( T, GTAB( IS) , BARBl(IS), BARB2(IS), BARB3(IS), 

DCP(IS)) 

52. Go to iu6 
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Since 


S = G + BARBl.Ln T + BARB2.T + BARB3.T 2 + DCP, ? - ~ Ln ^ 
c 1 1 I 1 loJP 

0 


P - P exp 
o 


G. + BARBljLn T + BARB2jT + BARB3^T 2 + DCPj y - S 


R/J 


The pressure P can be obtained, using entropy S and temperature T 
in the case of (IND1 = -1; T, S, — P) , where B is the input 
entropy. 


53. A = ( POMAX) ( Exp( ( GTAB( IS) + ( BARB1 ( IS) ) ( Ln( T) ) 
+ ( BARB2( IS) ) ( T) + ( BARB3( IS) ) ( T) 2 
+ ( DCP( IS) ) ( T) 3 / (3.0) - B)/ ROJ) ) 


54. Go to 106 

Variable specific heat is considered from 55, when IND1 = 2. 
Temperature T and specific heat C are obtained using given 
enthalpy H = (B). P 

(IND1 = 2; H — T, CP) 

Check whether the input enthalpy exceeds HTAB(IS) to search for 
temperature T e> 

55. If B ^ HTAB(IS), go to 59 
If B < HTAB(IS), go to 56 

56. IS = IS - 1 
Check the value of IS. 

57. If IS ^ 0, go to 42 
If IS > 0, go to 58 
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58. Go to 55 


59. If B < I1TAB(IS + l), go to (511 
If B a HTAB(IS + I), go to on 

00. IS = IS + 1 

01. If IS a 20, go to 42 
If IS < 20, go to 02 

02. Go to 59 

Check the value of IS. 

03. If IS — NOCTAB or IS s o, go to 42 

If IS < NOCTAB and IS > 0, go to 64 

Determine the table temperature and enthalpy corresponding to IS 
and IS + 1. 

64. TTP = TCTAB(IS) 

65. FP * HTAB(IS) 

66. TTPP * TCTAB(IS + l) 

67. FPP = HTAB(JS + 1) 

Determine the following quantity. 

68. TT0 = ((TTP) (FPP - B) - (TTPP)(FP - B))/(FPP - FP) 
Check the option indicator IND1. 

69. If IND1 > 2, go to 73 
If IND1 ^ 2, go to 70 
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Define DKLT (IND1 « 2; II - T, CP). 

70. DKLT = TTO - TCTAB(IS) 

Compute enthalpy. 

7 1 . FO - IITAB( IS) + ( C PTA B( IS) ) ( DE LT) + ( 0 . 5) ( PC P( IS) ) 

x (DEBT) 2 ! (CC PUS)) (DEBT) 3 / (3.0) 

+ (0.25) (DC P(IS)) (DEBT) 4 

72. Go to 74 

Compute entropy for IND1 = 0 (P, S -* T). 

73. FO = GAPF(TT0,GTAB(IS),BARB1(IS),BARB2(IS),BARB3(IS), 

DC P(IS)) 

Determine the temperature, 

74. TT1P = ( (TTO) (FPP - B) - (TTPP)(F0 - B))/(FPP - FO) 

75. TT1PP = ( (TTO) (FP - B) - (TTP)(F0 - B))/(FP - FO) 

Set switch indicator to minus one. 

76. N = -1 
Save temperature. 

77. TAU = TTO 
Save entropy. 

78. SF = FO 

Check whether the entropy or enthalpy calculated is nearly equal to 
the input one, using a tolerance of 10“ 7 . 

79. If I (SF - B) /(B) | ^ 10" 7 , go to 84 
If I (SF - B)/(B) | > 10" 7 , goto 80 
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Check whether the entropy o~ enthalpy computed is equal to or less 
than the input value. 

80. If SF S B, go to 86 
If SF > B, go to 81 

Set 

81. TTPP = TAU 

82. FPP = SF 
83« Go to 88 

Save the temperature in case convergence is achieved in step 79. 

84. T = TAU 

85. Go to 103 

Store temperature and entropy or enthalpy in case convergence is not 
obtained and the computed entropy or enthalpy is less than the input 
values. 

86. TTP = TAU 

87. FP = SF 

Check the switch indicator. 

88. If N < 0, go to 89 
If N = 0, go to 92 
If N > 0, go to 101 

Set the switch indicator to zero. 

89. N = 0 

Store the temperature. 


4 * 
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!)1. Go to 94 
Set N equal to one. 

92. N 1 

Store temperature TAU. 

92. TAU = TT1PP 
Chock temperature TAU, 

94. If TAU s TTP or TAU a TTPP, go to 88 
If TAU > TTP and TAU < TTPP, go to 95 

Check the option indicator IND1. 

95. If IND1 S 2, go to 98 
If IND1 > 2, go to 96 

Calculate entropy using temperature TAU, (IND1 = 3; P, S -*• T). 

96. SI*’ = GAPF( TAU, GTAB(IS) , BARBl(IS), BARB2(1S), 

BARB3( IS) , DC P(IS)) 

Repeat the iteration. 

97. Go to 79 

Define DELT to obtain enthalpy at temperature TAU 
(IND1 - 2; H — T, CP). 

98. DELT - TAU - TCTAB(IS) 

Enthalpy at temperature TAU; 

99. SI HTAB(IS) + ( CPTAB( IS) ) ( DELT) + ( 0. 5) ( BCP( IS) ) ( DELT) 

+ (CCP(IS)) (DELT) 3 / (3. 0) + ( 0. 25) ( DCP( IS)) (DELT) 4 



100. Go to 79 


Check whether convergence is obtained. 

101. If (FPP - FP)/(FP + FPP) > 0.0010, go to 68 
If (FPP - FP)/(FP + FPP) ^ 0.0010, go to 102 

Calculate the temperature. 

102. T = ( (TTP) (FPP - B) - (TTPP) ( FP - B))/(FPP - FP) 
Check the option indicator. 

103. If IND1 > 2, goto 106 
If IND1 ^ 2, go to 104 

Compute DELT, (IND1 =2, H — • T , CP). 

104. DELT = T - TCTAB(IS) 

Compute specific heat at temperature T. 

105. A = CPTAB(IS) + (BCP(IS))(DELT) + (CCP(lS))(DELT) 2 

+ (DCP(IS)) (DELT) 3 

Save the temperature. 

106. AA = T 

Save the specific heat or pressure. 

107. BB = A 

Check the option indicator. 

108. If IND1 = 2, return 

If IND1 * 2, go to 109 
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Save enthalpy or entropy depending on the option indicator IND1. 

109. CC = B 

110. Return 

Subroutine SEVAL nomenclature is given in Table 14. 


TABLE 14. SUBROUTINE SEVAL NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

A 

Saved value of 

pressure 

(JND1 = -1, 0,3) 

lbf/ft 2 

3, 16, 42, 50, 
53, 107 


or specific heat 
(IND1 = 1, 2) 

Btu/( lbm-° R) 

3, 9, 13, 105, 
107 

< 

Argument of 
entropy function 

Btu/(lbm-°R) 

1 

AA 

Saved value of 
temperature 

°R 

CALL, 2, 106 


Saved value of 
entropy 

(IND1 = .1, 0, 3) 

Btu/( lbm-° R) 

4, 18, 23, 42, 
53, 79, 80, 109 


Saved value of 
enthalpy 

(ft 2 /sec 2 ) or 

4, 6, 8, 14, 42, 
47, 48, 55, 59, 


(IND1 = 1, 2) 

Btu/ ( lbm- 6 R) 

68, 74, 75, 79, 
80, 109 

B1 

Argument of 
entropy function 

Btu/Ubm^R 2 ) 

1 

BARB1 

Polynomial 

equation for 

C = C . - BCP.T. 
P pi i i 

+ CCP^ 2 

- DCP.T? 

1 1 

Btu/(lbm-*R) 

/CSEVAL/, 17, 
22, 29, 51, 53, 
73, 96 

BARB2 

First derivative of 
BARB1 with negative 
sign 

Btu/ (lbm-°R 2 ) 

/CSEVAL/, 17, 
22, 29, 51, 53, 
73, 96 

BARB3 

Second derivative of 
BARB2 times - 1/4 

Btu/(lbm-°R 3 ) 

/CSEVAL/, 17, 
22, 29, 51, 53, 
73, 96 
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TABLE 14. SUBROUTINE SKVAL NOMENCLATURE (Continued) 



Description 


Specific heal in the 
case' of IND1 = 1 or 
2 

Pressure in the ease 
of IN 1)1 = -1, O ' or :i 


Coefficient in the 
C - T relationship 


Units 


Btu/(lbm-° R) 


lbf/ft 2 


( Btu/lbm) -° IP- 


determined by sub- 
routine BMFITS 


Argument of entropy 
function 


Enthalpy in the case 
of IND1 = 1 or 2 
Entropy in the case 
of IND1 = -1, 0 or 3 
Coefficient in the 
C - T relationship 


Btu/( lbm-° 


Reference 


CALL, 3, 


/CSEVAL/, 47, 
71, 99, 105 



CALL, 4 
CALL, 109 


/CSEVAL/, 47, 
71, 99, 105 



determined by sub- 
routine BMFITS 



CJG 

Specific heat at 
stagnation condition 
in work units 

ft 2 /( sec 2 -° R ) 

/CSEVAL/, 14 

CPO 

Specific heat at 
constant pressure at 
stagnation condition 

Btu/( lbm-° R) 

/CSEVAL/, 8, 
9, 13 

CPTAB 

Array of C values 
P 

corresponding to the 
values in the tem- 
perature table 

Btu/(lbm-°R) 

/CSEVAL/, 71, 
99, 105 

D1 

Argument of entropy 
function 

Btu/(lbm-°R 4 ) 

1 

DCP 

Coefficient in the 

C - T relationship 
P 

determined by sub- 
routine BMFITS 

Btu/(lbm-°R 4 ) 

/CSEVAL/, 17, 
22, 29, 47, 51, 
53, 71, 73, 96, 
99, 105 
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TABLE 14. SUBROUTINE SEVA L NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

BELT 

Temperature 

difference 


40, 47, 70, 71, 
08, 00, 104, 105 

FO 

Saved value of 
enthalpy 
Saved value of 
entropy 

Ulu7ibin 

IHu/( lbm-° It) 

71, 74, 75, 78 
73, 74, 75, 78 

FJG 

Conversion factor 
between thermal and 
work units multiplied 
by acceleration of 
gravity used as a 
proportionality 
constant 

( fL 2 -lbmj / 
(Btu-soc 2 ) 

/CSKVAL/, 0, 
48 

FP 

Saved value of 
table enthalpy 

ft 2 /sec 2 

65, 68, 75, 87, 
101, 102 

FPP 

Saved value of table 
enthalpy 

ft 2 / sec 2 

67, 68, 74, 82, 
101, 102 

G 

Argument of entropy 
function 


1 

GAPF 

Entropy function 

Btu/( lbm-° R) 

1, 17, 22, 29, 
51, 73, 96 

GTAB 

Array of terms used 
in entropy function 

^Btu/( lbm-° R) 

/CSEVAL/, 17, 
22, 29, 51, 53, 
73, 96 

HTAB 

Array of enthalpy 
values 

Btu/lbm 

/CSEVAL/, 47, 
55, 59, 65, 67, 
71, 99 

IND1 

Option indicator 


CALL, 5, 11, 
42, 45, 69, 95, 
103, 108 

ICTAB 

Indicator = 0, 

constant specific 

heat calculation 

= (3 ^ ICTAB s 20), 

dimension values in 

C versus T tables 
P 


/INPUT/, 7, 12 
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TABLE 14. SUBROUTINE SEVAL NOMENCLATURE (Continued) 



Symbol 

Description 

Units 

Reference 


IS 

Tabulated position 

_ 

/CSEVAL/, 17, 



in flow direction for 


19, 20, 22, 24, 



temperature table 


25, 27-30, 33- 





35, 37-39, 41, 





42, 44, 46, 47, 





51, 53, 55-57, 

I 




59-61, 63-67, 

l 




70, 71, 73, 96, 

1 




98, 99, 104, 105 

| 

N 

Switch indicator 

— 

76, 88, 89, 92 

■I 

NOCTAB 

Saved value of 

— 

/CSEVAL/, 27, 

[ 


ICTAB 


44, 63 

1 

POMAX 

Saved value of 

lbf/ft 2 "™ 

/CSEVAL/, 16, 

\ 


stagnation pressure 


50, 53 


PR 

R ' ln (t) 

Btu/( lbm-° R) 

1, 16, 50 


ROJ 

Specific gas constant 

Btu/(lbm-° R) 

/CSEVAL/, 16, 

5 


in thermal units, 


50, 53 

* 4 


RBAR divided by FJ 



* 

SF 

Saved value of 

Btu/lbm 

78-80, 82, 87, 



enthalpy or 


99 

j 


Save value of 

Btu/(lbm-° R) 

78-80, 96 

. 


entropy 



1 

STAB 

Computed value of 

Btu/(lbm-° R) 

17, 18, 22, 23 

V 


entropy 




T 

Temperature 

°R 

1, 2, 8, 14, 33, 



1 


37, 42, 46, 51, 

?: 




53, 84, 102, 





104, 106 


TAU 

Saved value of 

°R 

77, 81, 84, 86, 

- 


temperature 


90, 93, 94, 96, 

* 




98 

-- 

TCTAB 

ICTAB values of 

°R 

/CSEVAL/, 17, 



temperature in 


22, 28-30, 33, 



increasing order 


37, 46, 64, 66, 



defining C - T 


70, 98, 104 



table p 
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TABLE 14. SUBROUTINE SEVAL NOMENCLATURE (Concluded) 


Symbol 


Description 


Units 


Reference 


TTO 


TT1P 


TT1PP 


Assumed te mpera ture 
for iteration on the 
tc mpera ture 


R 


A ssumed tempera tu re 
for iteration on the 
temperature 


°R 


Assumed temperature 
for iteration on the 
temperature 


R 


(58, 70, 73-75, 
77 


74, 90 


75, 93 


TTP 

TTPP 


Saved value of 
TCTAB (IS) 


'R 


64, 68, 75, 86, 
94, 102 


Saved value of 
TCTAB (IS + 1) 


R 


66, 68, 74, 81, 
94, 102 
















SUBROUTINE START 


Subroutine START locates the sonic point (M ~ 1) from the input Mach 
number versus x table and then computes initial values for momentum and 
energy thicknesses. 

The table of Mach numbers is first scanned to locate an exact input 
value of M = 1. If none is found, the sonic point location is determined from 
the following iteration procedure. 

1. Locate the i^ interval, in the table, which contains M = 1, 

and assume x = x. . 

g i 

2. Using subroutine XNTEltP. evaluate M at x . 

g g 

3. Make a tolerance comparison. 

If |M - 1.0 |s TOLZME, convergence is satisfactory, 
g 

and appropriate flow quantities are avaluated at x . If 

g 

outside the tolerance, a secant method is used to make 
another assumption for x^ and steps 2 to 3 are repeated 

up to a maximum of 50 iterations. 

Using the flow quantities just computed, subroutine INTZET is used 
to evaluate integrals I t and I 2 for C, = 1 . The shape factor 6*/0 is then 
evaluated from equation (1251 11] . An iteration for the skin friction coefficient 
is performed (similar to that in subroutine BARPRO ) , and a value of 

momentum thickness e is computed. These values are assumed to be the 

initial values for 6 and <p at x . 

g 

COMMON BLOCKS 

COMMON blocks COFIIF, CSEVAL, INPUT, INTER, 

LOOKUP, OUTPUT, SAVED, and TABLES are used. _ 

TBL SUBROUTINES 

Subroutine BARCON calls START . 

START calls CFEVAL, GETPT, INTZET, QUITS, 

SEVAL, and XNTERP . 
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FORTRAN SYSTEM ROUTINES 

FORTRAN library routine SQltT is used. 

Built-in FORTRAN function ABS is used. 

CALLING SEQUENCE 

The subroutine calling sequence is: 

CALL START 
SOLUTION METHOD 

Set table position indicators for a first time 
entry indication. 

1. IMX = _l 

2. IYS = -1 

3. ITWX = -1 

Initialize subscript counter to zero. 

4. I = 0 

Increment subscript counter. 

5. 1=1+1 

Check whether Mach number table includes M = 1.0 . 

6. If ZMTAB(I) < 1.0 , goto 7 
If ZMTAB(I) = 1.0 , go to 10 
If ZMTAB(I) > 1.0 , goto 13 

Check whether subscript counter is equal to the table dimension. 
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7. If I < IXTAB , go to 5 


If I IXTAB , go to 8 

The Mach number table does not contain M - 1.0 . 

8. WRITE error message about Mach number table 
Print out all the COMMON blocks and go to next case. 

0. CALL QUITS 

Obtain axial distance at which M = 1. 0 . 

10. X - XITAB(I) 

Set subscript counter for the Mach number that just exceeds one. 

11. IBEG =1+1 

12. Go to 29 

Check whether the first Mach number table value is greater than one. 

13. If I ^ 1 , go to 8 
If I > 1 , go to 14 

Save value of axial distance at which Mach number becomes greater 
than one. 

14. XG = XITAB(I) 

Obtain the lowest number table value greater than one. 

15. ZME = ZMTAB(I) 

Approximate the axial distance at which M = 1. 0 . 

16. X = XITAB(I-l) + (XITAB(I) - XITAB(I - 1))/ 

(ZMTAB(I) - Z MTAB (I - l) ) (l. 0 - ZMTAB(l - l)) 

Initialize Mach number iteration counter to zero. 



Set subscript counter for which Mach number exceeds one. 

18. IBEG - I 

Increment Mach number iteration counter, 

19. J = J + 1 

Save previous value of axial distance at which M > 1.0, 

20. XO = XG 

Save previous value of Mach number obtained from table lookup. 

21. ZMO = ZME 

Save present value of axial distance for which M > 1.0. 

22. XG » X 

Determine value of Mach number for approximated value of axial 
distance. 

23. CALL XNTERP(X, ZME . ZMEP . IMX, XITAB, ZMTAB, KTAB 

CMX, IMX) 

Check whether the Manh number falls within the desired tolerance 
band. 

24. If | ZME - 1.0| TOLCFA , goto 29 
If | ZME - 1.0| > TOLCFA, go to 25 

Compute a new value of the Mach number. 

25. ZMX = (XG - XO)/ (ZME - ZMO) 

Compute a new value of the axial distance. 




26. X = XO + (1.0 - ZMC)) (ZMX ) 

Check whether there have been more than 60 iterations to determine 
the Mach number equal to one. 

27. It J GO, goto 19 
It J > GO , go to 28 

There have been more than 50 iterations on the Mach number. 

28. WRITE error message about Mach number iteration failure 
Determine the value of Mach number at final value of axial distance. 

29. CALL XNTERP (X, ZME , ZMEP, IMX, XITAB, ZMTAB, IXTAB, 

CMX, IMX) 

Determine the nozzle radius at this value of axial distance. 

30. CALL XNTERP (X, YR, YRP, IYS, XITAB, YITAB, IXTAB, 

CYX, IMX) 

Compute the pressure and temperature for the computed Mach number. 

31. CALL GETPT (ZME, PSE, TE) 

Compute the enthalpy and specific heat at this temperature. 

32. CALL SEVAL (l, TE, CPE, HE) 

Compute the ratio of specific heats. 

33. GAME = (CPE) /'(CPE - (RBAR)/(FJ)) 

Compute the difference between stagnation and static enthalpy. 

34. HB - HO - HE 
Compute the fluid velocity. 
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35. UE = SQRT( (2. 0) (HB) ) 

Compute the fluid density. 

36. RHSE = (PSE)/(TE)/(RBAR) 

Compute viscosity at the static temperature. 

37. ZMU = (ZMU0)((TE)/(T0)) ZMVIS 

Compute the adiabatic wall enthalpy. 

1 / 

38. HAW = HE + ((PRANDT) 3 ) (HB) 

Compute the adiabatic wall temperature and specific heat. 

39. CALL SEVAL(2, TAW . CPAW . HAW) 

Check for temperature option to be used. 

40. If ITWTAB < 0 , go to 41 

If ITWTAB = 0 , go to 44 

If ITWTAB > 0 , go to 47 

Adiabatic wall temperature option is used; set wall enthalpy and wall 
temperature to adiabatic values. 

41. HW = HAW 

42. TW - TAW 

43. Go to 49 

Constant wall temperature option is used; set wall enthalpy to previously 
computed value. 

44. HW = TWTAB(2) 

Set wall temperature to input value. 
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45. TW ... TWTAJi(l) 


46. Go to 49 

Tabular wall option is used; determine the value of wall temperature. 

47. CALL XNTEHP(X, TW, TW1>, 1TWX, XITAB, TWTAU, IX TAB, 

C'J'WX, MX) 

Compute enthalpy and specific heat at the wall temperature. 

48. CALL SEVA L(l, TW, CPW , 11W) 

Obtain flow properties to evaluate integrals Ij and I 2 . 

49. AFINT - HW 

50. BFINT = HO - HW 

51. CFINT = - HB 

52. TFINT - TE 

53. MMINT = MZETA 

Set indicator to evaluate first integral Ij . 

54. IFINT = 1 
Evaluate first integral Ij . 

55. CALL INTZET( 0.0, 1.0, ZI1) 

Set indicator to evaluate second integral I 2 . 

56. IFINT = 2 

Evaluate second integral I 2 . 

57. CALL INTZETi 0.0, 1.0, ZI2) 

Compute the boundary layer shape factor. 

i 

* 

i 
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58. DELSOT - ((l.O)/(ZM/KTA) - ZI2)/(ZH) 

Check for type of flow to be considered. 

59. If KPSZ = 0.0, go to 02 
If KPSZ + 0.0, go to 00 

Axisymmctric flow is considered; set intermediate term to ratio of 
nozzle radius slope to the nozzle radius. 

00. ERASES = (YRP)/(YR) 

01. Go to 03 

Two-dimensional flow is considered; set intermediate term to one. 

62. ERASE5 =1.0 
Compute intermediate term. 

63. ERASE4 = (l.O + DELSOT) /( 1.0 + (GAME - 1. 0) /(2. 0) ) (ZMEP) 

+ ERASE 5 

Check on value of intermediate term. 

64. If ERASE4 * 0.0, go to 66 
If ERASE4 = 0.0, go to 65 

The intermediate term is in error. 

65. WRITE error message and possible remedy 

Compute ratio of difference between stagnation and adiabatic wall 
temperature to the adiabatic wall temperature. 

66. ERASE 1 = (17. 20) (TO - TAW)/ (TAW) 

Compute ratio of difference between static and stagnation temperature 
to adiabatic wall temperature. 


i 
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67. ERASE2 = (305.0)(TE - TO) /(TAW) 

Compute term for calculating the momentum thickness. 

68. CTHET ~ ( 0. 50) (SQRT( 1. 0 + YRP 2 ) ) /( ERASE4) 

Compute skin friction coefficient intermediate term. 

59. CRT2 - ((TAW)/(TE))^ 1,0 " ZMVIS )(RHSE)(UE)/ 

(<ZMU) (CTHET) 

Compute ratio of adiabatic wall temperature to the static temperature. 

70. ERASE3 = (TAW)/(TE) 

Set initial assumption for the skin friction coefficient. 

71. CFA = 0.0010 

Initialize skin friction coefficient iteration counter to zero. 

72. JE = 0 

Increment skin friction coefficient iteration counter. 

73. JE = JE + 1 

Save previous skin friction coefficient quantity 

74. CFG = CFA 

Compute momentum thickness based on computed skin friction 
coefficient. 

75. THETA = (CFG) (CTHET) 

Compute term for obtaining skin friction coefficient based on Coles’ 
relationship. 

76. CR = (CRT2) (THETA) 2 
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Compute skin friction coefficient based on Coles' relationship. 

77. CFB = CFEVAL(CR) 

Compute intermediate term. 

78. TTAW -- 1.0+ ( ERASED (SQIIT(( CFB)/ (2.0) ) 

+ ( ERASE2) (CFB)/ ( 2, 0) 

Check value of intermediate term. 

79. If TTAW > 0.0, go to 86 
If TTAW =£ 0.0, go to 80 

Compute new value of the skin friction coefficient. 

80. CFA = (3.0) (CFG) 

Save computed value of skin friction coefficient. 

81. Z4 = CFA 

Save approximation value of skin friction coefficient. 

82. Z2 = CFG 

Check whether there haye been over 50 integrations to determine the 
skin friction coefficient. 

83. If JE 50, go to 73 
If JE > 50, go to 84 

There have been more than 50 iterations on the skin friction 
coefficient. 

84. WRITE error message about skin friction coefficient iteration 
failure. 

85. Go to 96 




f 


Compute new value of the skin friction coefficient. 

86. CFA = (CFB) /((ERASES) (TTAW) ZMVIS ) 

Check whether the skin friction coefficient is within the tolerance. 

87. If ! (CFA - CFG) /(CFA + CFG) | ss TOLCFA, go to 96 
If | (CFA - CFG) /(CFA + CFG) | > TOLCFA, go to 88 

Check whether the iteration counter is less than two. 

88. If JE < 2, go to 81 
If JE ^ 2, go to 89 

Save previous computed value of skin friction coefficient. 

89. Z3 = Z4 

Save previous assumed value of skin friction coefficient. 

90. Z1 = Z2 

Save new computed value of skin friction coefficient. 

91. Z4 = CFA 

Save new approximation of skin friction coefficient. 

92. Z2 = CFG 

Com pute skin fri ction coefficient convergence term. 

93. ZS5 = (Z4 - Z3)/(Z2 - Zl) 

Compute new value of the skin friction coefficient. 

94. CFA = (Z4 - (ZS5)(Z2))/(1.0 - ZS5) 

Loop back to continue iteration on skin friction coefficient. 

95. Go to 83 
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Compute initial value of momentum thickness. 

96. THETAI = (CFA)(CTHET) 

Set initial value of momentum thickness equal to energy thickness. 

97. PHII = THETAI 

Save final computed value of skin friction coefficient. 

98. CFAGT = CFA 
Set shape factor to one. 

99. ZETA =1.0 

Write out the initial values of momentum and energy thicknesses at the 
throat. 

100. WRITE X, YR, THETAI, PHII 

101. Return 

Subroutine START nomenclature is given in Table 15. 


TABLE 15. SUBROUTINE START NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

AFINT 

Wall enthalpy for 
integral evaluation 

ft 2 /sec 2 

/COFIIF/, 49 

BFINT 

Difference between 
stagnation enthalpy 
and wall enthalpy 

ft 2 /sec 2 

/COFHF/, 50 

CFA 

Computed value of 
skin friction 
coefficient 


71, 74, 80, 81, 
86, 87, 91, 94, 
96, 98 

CFAGT 

Saved value of 
computed skin 
friction coefficient 

” 

/INTER/, 98 

CFB 

Value of skin friction 

coefficient 

from CFEVAL 


77, 78, 86 










TABLK 15. SUBROUTINE START NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

CFG 

Approximation of 
akin friction 
coefficient 

— 

74, 75, 80, 82, 
87, 92 

CFINT 

Difference between 
static enthalpy 
and stagnation 
entha lpy 

ft 2 /sec 2 

/COFIIF/, 51 

CMX 

Array of parabola 
coefficients for 
Mach number table 
(ZMTAB) 


/LOOKUP/, 23, 
29 

CP AW 

Specific heat at 
adiabatic wall 
conditions 


39 

CPE 

Specific heat at 
static temperature 

Btu/( lbm-° R) 

32, 33 

CPW 

Specific heat at 
wall temperature 


48 

CR 

Term for computing 
skin friction 
coefficient in 
CFEVAL 

lbm/(ft-sec) 

76, 77 

CRT2 

Term for 
computing CR. 

lbm/(ft 3 -sec) 

69, 76 

CTHET 

Term for computing 
momentum thickness 

ft 

68, 75, 96 

CTWX 

Array of parabola 
coefficients for the 
wall temperature 
table (TWTAB) 


/LOOKUP/, 47 

CYX 

Array of parabola 
coefficients for the 
nozzle radius table 
(YITAB) 


/LOOKUP/, 30 


Boundary layer shape 
factor 6*/0 


/OUTPUT/, 58, 

63 , 


Flow geometry 
indicator 




189 
































TABLE 15. SUBROUTINE START NOMENCLATURE (Continued) 


Symbol 

Description 

ERASE 1 

Ratio of stagnation 
and adiabatic wall 
temperature 

ERASE 2 

Ratio of static, 
stagnation, and 
adiabatic wall 
temperatures 

ERASE 3 

Ratio of adiabatic wall 
temperature to the 
static temperature 

ERASE4 

Intermediate term 
for computing 
momentum thickness 

ERASE 5 

Intermediate term 
for computing 
momentum thicliness 

FJ 

Thermal-to-wo:’k 
unit conversion 
factor 

GAME 

Ratio of specific 
heats 

HO 

Stagnation enthalpy 

HAW 

Adiabatic wall 
enthalpy 

HB 

Difference between 
stagnation and 
static enthalpy 

HE 

Static enthalpy 

HW 

Wall enthalpy 

I 

Array subscript 
counter 

' IBEG 

Subscript counter at 
which the Mach 
number table exceeds 
one 


Units 


ft 2 /sec 2 


ft 2 /sec 


ftVsec 


Reference 
66, 78 

67, 78 
70, 86 


63, 64, 68 


60, 62, 63 


/INPUT/, 33 


33, 63 


/CSEVAL/, 34, 
50 


38, 39, 41 


34, 35, 38, 51 


/INTER/, 32, 
34, 38 


/INTER/, 41, 
44, 48-50 


4-7, 10, 11, 13- 
16, 18 


/INTER/, 11, 18 




























TABLE 15. SUBROUTINE START NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

IFINT 

Integral evaluation 
indicator 

* 

/COFIIF/, 54, 
56 

IMX 

Position or start 
indicator for Mach 
number table 
(ZMTAB) 


/LOOKUP/, 1, 
23, 29, 30, 47 

ITWTAB 

Wall temperature 
indicator 

— 

/INPUT/, 40 

ITWX 

Position or start 
indicator for wall 
temperature table 

(twtab) 

t; 

/LOOKUP/, 3, 
47 

IXTAB 

Nvmber of points in 
x, y, and M tables 


/INPUT/, 7, ” 

23, 29, 30, 47 

IYS 

Position or start 
indicator for nozzle 
radius table (YITAB) 


2, 30 

J 

Mach number 
iteration counter 

— 

17, 19, 27 

JE 

Skin friction 
coefficient iteration 
counter 


72, 73, 83, 88 

MMINT 

Exponent for integral 
evaluation 

— 

/COFIIF/, 53 

MZETA 

Exponent of velocity 
profile 

— 

/INPUT/, 53 

PHII 

Initial value of 
energy thickness 

ft 

/INPUT/, 97, 
100 

PRANDT 

Prandtl number 

— 

/INPUT/, 38 " 

PSE 

Static pressure at 
M = 1.0 

lbf/ft' 1 

31, 36 

RBAR 

Specific gas constant 

( ft-lbf) / ( lbm-° R) 

/INPUT/, 33, 36 

RHSE 

Fluid density at 
M = 1.0 

lbm/fr 

36, 69 

TO 

Stagnation 

temperature 

0 R 

“7lNPUT/, 37, 
66, 67 

TAW 

Adiabatic wall 
temperature 

°R 

39, 42, 66, 67, 
69, 70 
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TABLE 15. SUBROUTINE START NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

TE 

Static temperature 
at M = 1.0 

°R 

/OUTPUT/, 31, 
32, 36, 37, 52, 
67, 69, 70 

TFINT 

Static temperature 
for integral, 
evaluation 

0 R 

/COFIIF/, 52 

THETA 

Boundary layer 
momentum thickness 

ft 

/OUTPUT/, 75, 
76 

THETAI 

Initial value of 
momentum thickness 

ft 

/INPUT/, 96, 
97, 100 

TOLCFA 

Tolerance in skin 
friction coefficient 
iteration 


/INPUT/, 24, 
87 

TTAW 

Intermediate term 
in skin friction 
coefficient iteration 


78, 79, 86 

TW 

Wall temperature 

° R 

/OUTPUT/, 42, 
45, 47, 48 

TWP 

Derivative of wall 
temperature 

0 R/ft 

47 

TWTAB 

Array of wall 
temperatures at 
contour points 

°R 

/TABLES/, 47 

UE 

Velocity of fluid 
at M = 1.0 

ft/sec 

/OUTPUT/, 35, 
69 

X 

Axial distance at 
which M = 1.0 

ft 

/OUTPUT/, 10, 
16, 22, 23, 26, 
29, 30, 47, 100 

XG 

Approximation of 
axial distance at 
which M = 1.0 

ft 

14, 20, 22, 25 

XITAB 

Array of axial 
distance values 

ft 

/TABLES/, 10, 
14, 16, 23, 29, 
30, 47 

XO 

Previous assumption 
of axial distance at 
which M = 1.0 

ft 

20, 25, 26 

YITAB 

Array of nozzle 
radius values 

ft 

/TABLES/, 30 
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TABLE 16. SUBROUTINE START NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

YR 

Nozzle radius or 
contour height 

ft 

/OUTPUT/, 30, 
60, 100 

YRP 

Derivative of nozzle 
radius 



30, 60, 68 

Z1 

Saved value of skin 
friction coefficient 
assumption 


/OUTPUT/, 90, 
93 

Z2 

Saved value of skin 
friction coefficient 
assumption 


/OUTPUT/, 82, 
90, 92-94 

Z3 

Saved value of 
computed skin 
friction coefficient 


/OUTPUT/, 89, 
93 

Z4 

Saved value of 
computed skin 
friction coefficient 


/OUTPUT/, 81, 
89, 91, 93, 94 

ZETA 

Shape factor 
(A/5 ) 1/n 


/OUTPUT/, 99 

zii 

Value of integral I t 

- 

/SAVED/, 55, 
58 

ZI2 

Value of integral I 2 

— 

/SAVED/, 57, 58 

ZME 

Mach number 

1 

/OUTPUT/, 15, 
21, 23-25, 29, 31 

ZMEP 

Derivative of Mach 
number 

ft -1 

23, 29, 63 

ZMO 

Saved value of Mach 
number for iteration 
on M* 1.0 


21, 25, 26 

fclVITAB ' 

Array of Mach 
numbers associated 
with contour points 


/TABLES/, 6, 
15, 16, 23, 29 

ZMU 

Fluid viscosity at 
M = 1.0 

lbm/(ft-sec) 

37, 69 

ZMUO 

Viscosity at 

stagnation 

temperature 

lbm/(ft-sec) 

/INPUT/, 37 

ZMVIS 

Exponent in viscosity 
temperature law 


/INPUT/, 69, 86 





Symbol 


Description 


Reference 


ZMX Computed value of 

Mach number in 

iteration on M = 1. 0 

ZMZETA Real value of 

velocity profile 
exponent 

Ratio of difference 
between computed 
values of skin 
friction coefficient 
and approximations 


25, 26 


/INTER/, 58 


93, 94 
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SUBROUTINE XNTERP 


Subroutine XNTERP determines the dependent variable y and its 
first derivative y' for a dependent variable x based upon a functional 
relationship y = f(x) represented in tabular form. A local quintic spline 
interpolation is used in the solution process as described below. 

The independent variable table is searched until the interval containing 

x is located. The associated interval table values x, and x and two 

1 1 + 1 

additional ones x, , and x. , _ , one on either side of the specified interval, 

1+2 

are then selected for the quintic spline fit. Using either the first three 
x. x , x., x. + jor the last three^x., x. + 1> x. +2 jtable values from the ones 

previously specified, two parabolas and their corresponding coefficients are 
determined satisfying the following relationship. 

Y l = Ci + C 2 (x - x ul ) + C 3 (x - x ul ) 2 , 

Y r = C 4 + C 5 ( X " X i) + C 6 (X - X.) 2 , 

= C 2 + 2C 3 (x - x ul ) , 

7^ = C 5 + 2C 6 (x - x.) . 


RIGHT PARABOLA 



A cubic weighting function a and its derivative a * 
a (x) = 3U 2 - 2U 3 , 

a ' (x) = (6U - 6U 2 ) • U'(x) , 
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provide a scheme to solve for the dependent variable y and its derivative y' 
according to the relationship below. 


y (x) = [ 1 - a (x)] [Cj + C 2 (x - x.^) + C 3 (x - x ui ) 2 ] 

+ a (x) [C 4 + C 5 (x - x^) + C 6 (x - x.) 2 ] 

■ tl - « (x) ] y L + a (x) f R ; 

y’(x) = a* (x)(7 r - y L ) + [1 - a (x) 1 y^+ a (x) . 

COMMON BLOCKS 

No COMMON blocks are used. 

TBL SUBROUTINES 

Subroutines BARCON . BARPRO. and START call XNTERP. 
XNTERP calls subroutine QUITS. 

FORTRAN SYSTEM ROUTINES 

No FORTRAN library routines or built-in FORTRAN functions 
are used. 

CALLING SEQUENCE 

The subroutine calling sequence is 

CALL XNTERP (X, X, X£, DON, XAR, YAR, 1A.% CAR, IPOS) 



X = input value of the independent variable 
Y — output value of the corresponding dependent variable 
YP = derivative of the dependent variable 
IXIN = table position or start indicator 
XAR = input array of independent variables 
YAR = input array of dependent variables 
jar = dimension of the independent and dependent arrays 
CAR = input or output array of parabola coefficients 
IPOS = array position indicator 
SOLUTION METHOD 

Set interval subscript counter Limit. 

1. IXO = IXIN 

Set working maximum length of input arrays. 

2. EX MAX = EAR - 1 

Set internal subscript counter to array position indicator. 

3. IX = IPOS 

Set array of parabola coefficients to those used previously, 

4. Do 5, I = 1, 6 

5. C ( I) = CARl I) 



Check whether this is a first time entry. 

6. If IXO > 0, go to 10 
If IXO 0, go to 7 

This is a first time entry} set the indicator to one, 

7. IF1RST = 1 

Set subscript counter limit to IAR + 1. 

8. IXO = IX MAX + 2 
Set array subscript counter to one. 

9. DC = 1 

Check whether subscript counter is zero or negative. 

10. If IX < 0, go to 7 

If DC > 0, go to 11 9 

Check whether the input independent value is greater than or equal to | 

table value. 

11. If X a XAR(IX) , go to 16 
If X < XAR(DC), go to 12 

Decrement the array subscript counter. 

12. IX = DC - 1 

Check whether the subscript counter is still greater than zero. 

13. If IX > 0, go to 11 
If IX s 0, go to 14 



! 7.1 ? }< 



i 


t 


The input value is out of range of the table; print out values. 


14. WRITE X, XAR( 1) , XAR( IXMAX + 1) , YAR( 1) , 
YARdXMAX + 1) 


Print out the COMMON blocks, and go to the next case. 

15. CALL QUITS 


Check whether the input value is less than or equal to the next table 
value. 


16. If X < XAR(IX + 1) , go to 19 
If X > XAR(IX + 1) , go to 17 

Increment the array subscript counter. 

17. K = K + 1 


Check whether the subscript counter has exceeded the maximum value. 

18. If DC < IXMAX, go to 16 
If IX > IXMAX, go to 14 


The points bracketing x have been found; obtain the four surrounding 
points for x and y . 


19. Do 22, 1=1,4 


Set subscript counter to obtain x^ , x^ , x^, and x^ + ^ • 


20. II = IX - 2 + I 


Obtain the four surrounding points from the x table and the corres- 
ponding points from the y table. 


21. XI(I) = XAR(Il) 

22. YI(I) = YAR(Il) 
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A 

p 


Compute difference between x and x^. 


23. DX2 = X - XI( 2) 


Compute difference between 


i+1 


and x 


24. DX32 = XI( 3) - XI(2) 


Check subscript counter against subscript limit, 

25. If IX < IXO, go to 39 


If IX = IXO, go to 26 
If IX > IXO, go to 62 

The subscript limit has been reached; set subscript limit indicator. 

26. IXOGO = 0 


Check whether subscript counter is greater than one. 

27. If IX > 1, go to 30 
If IX s l, go to 28 

Subscript counter is less than or equal to one; set indicator. 

28. IGO = -1 

29. Go to 74 

Check whether subscript counter is less than the maximum. 

30. If IX < IX MAX, go to 37 
If EX ^ IX MAX, go to 31 

Check for first time entry into the dependent array. 

31. If IFIRST = 0, go to 35 
If IFIRST * 0, go to 32 


This is a first time entry; reset entry indicator, and set loop indicator. 

32. IFIRST = 0 

33. IGO =■- 1 

34. Go to 54 

This is not a first-time entry; set the loop indicator. 

35. IGO = 1 

36. Go to 70 

Subscript counter is less than the maximum; set loop indicator. 

37. IGO = 0 

38. Go to 70 

Subscript counter is below the limit; set subscript limit indicator. 

39. IXOGO = -1 

Check whether the subscript counter is one less than the limit. 

40. If IX < IXO -1, go to 45 
If IX > IXO - 1, go to 41 

Set left parabola coefficients to the right parabola coefficients. 

41. C(4) = C(l) 

42. C( 5) = C( 2) 

43. C(6) = C( 3) 

44. Go to 52 

Compute coefficients for the left parabola. 

45. C(4) = YI(2) 
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Compute difference between x and x. and y and y, . 

1+ — i i+ 1 i 

46. DX42 = XI(4) - XI(2) 

47. DY32 - YI( 3) - YI(2) 

48. DYOX32 - ( DY32)/ ( DX32) 

49. C(6) = ( DYOX32 - (YI(4). - „YIi2) )/ ( DX42) )/ (XI( 3) - XI(4)) 

50. C( 5) - DYOX32 - (C(6))(DX32) 

Check subscript limit indicator. 

51. If IXOGO > 0, go to 70 
If IXOGO s 0, go to 52 

Check whether subscript counter is greater than one. 

52. If IX s l, go to 28 
If IX > 1 , go to 53 

Subscript is less than the limit and greater than one; set indicator. 

53. IGO = 0 

Compute coefficients for the left parabola. 

54. C(l) = YI(1) 

Compute difference between x. and x^ . 

55. DX21 = XI( 2) - XI(1) 

Compute difference between x. and x 

l+l i-1 

i-1* 


56. DX31 = XI( 3) - XI(1) 

Compute difference between y^ and y 

57. DY2.1 - YI( 2) - YI( 1) 


58. DYOX21 = (DY2l)/(DX2l) 

59. C(3) = (DYOX21 - (Yl(3) - Yl( l))/(DX3l))/(XI(2) - XT(3)) 

60. C(2) = DYO X21 - (c(.3))(DX2l) 

Check subscript limit indicator. 

61. If IXOGO £ 0, go to 70 
If IXOGO > 0, go to 67 

The subscript counter is above the limit; set subscript limit indicator. 

62. IXOGO = 1 

Check whether the subscript counter is more than one above the limit. 

63. If IX > IXO + 1, go to 54 
If IX ^ IXO + 1, go to 64 

Set coefficients of left parabola to those of right parabola. 

64. Ctl) = C( 4) 

65. C(2) = C(5) 

66. C(3) = C(6) 

Check whether subscript counter is less than the maximum. 

67. If IX 2= IX MAX, go to 35 
If IX < EXMAX, go to 68 

Subscript counter is less than the maximum set loop indicator. > 

68. IGO = 0 

69. Go to 45 

f 
V 

.1 

l 
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Compute difference between x and x. . 

i -1 

70. DX1 = X - XI(1) 

Compute dependent var iable value from the left parabola. 

71. YB1 ■ ((C(3))(DX1) + C(2))(DX1) + C(l) 

Compute derivative of dependent variable. 

72. YPB1 = (C(3))(DXl)/(0.50) + C(2) 

Check loop control indicator. 

73. If IGO > 0, go to 89 
If IGO 0, go to 74 

Compute dependent variable value from the right parabola. 

74. YB2 = ((C(6))(DX2) + C(5))(DX2) + C(4) 

Compute derivative of dependent variable. 

75. YPB2 = (C(8))(DX2)/(0.50) + C(5) 

Check loop control indicator. 

76. If IGO < 0, go to 92 
If IGO a 0, go to 77 

Compute U (x) , U 2 , and U 3 for the cubic weighting function. 

77. U1 = (DX2)/(DX32) 

78. U2 - (U1)(U1) 

79. U3 = (U1)(U2) 

Compute the cubic weighting function a . 

80. Al = (3.0)(U2) - (2.0)(U3) 
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Compute the derivative of the cubic weighting function. 

81. A1P = (6.0)(U1 - U2)/(DX32) 

Compute the interpolated value of the dependent variable, 

82. Y « (1.0 - Al) ( YBl) + (A1)(YB2) 

Compute the derivative of the dependent variable, 

83. YP = (1.0 - Al) ( YPBl) - ( A1P) (YBl - YB2) + (Al)(YPB2) 

Save value of the array subscript counter for the next entry into the 
dependent table. 

84. DON = IX 

Check whether the subscript counter limit has been reached. 

85. If IXOGO = 0, return 
If IXOGO * 0, go to 86 

Save the values of the parabola coefficients for the next entry into the 
same tables. 

86. Do 87, I = 1, 6 

87. CAR(I) = C(I) 

Return to the calling subroutine. 

88. Return 

Set the value of the dependent variable to that computed from the left 

parabola. * 

89. Y = YBl 

Set the derivative to that computed from the left parabola, 

90. YP = YPBl 
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Save array position indicator and parabola coefficients. 

91. Go to 84 

Set the value of the dependent variable to that computed from the right 
parabola. 

92. Y = YB2 

Set the derivative to that computed from the right parabola. 

93. YP = YPB2 

Save array position indicator and parabola coefficients. 

94. Go to 84 

Table 16 gives the subroutine XNTERP nomenclature. 


TABLE 16. SUBROUTINE XNTERP NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

A1 

Cubic weighting function a (x) 


80, 82, 83 

A1P 

Derivative of the cubic weighting 
function a 1 (x) 


81, 83 

C 

Internal array of coefficients for 
right and left parabolas 


DIM, 5, 41-43, 
45, 49, 50, 54, 
59, 60, 64-66, 
71, 72, 74, 75, 
87 

CAR 

Input or output array of 
parabola coefficients for 
input tables 


CALL, DIM, 5, 
87 

DX1 

Difference between value of 
independent variable x and 
the second table point less than 
X (x . x ul ) 


70, 71 

DX2 

Difference between value of 
independent variable x and the 
first table point less than or 
equal to x (x - x ) 


23, 72, 74, 75, 
77 


206 












TABLE 16. SUBROUTINE XNTERP NOMENCLATURE (Continued) 




Symbol 

Description 

Units 

Reference 

DX21 

Difference between first 1 ble 
point less than x and the second 
table point less than x (x^ - x^ 


55, 58, 60 

DX31 

Difference between first table 
point greater than x and second 
table point less than 

x <Vi - x i.i> 


56, 59 

DX32 

Difference between first table 
point greater than x and first 
table point less than 
x (x i+1 - X,) 


24, 48, 50, 77, 
81 

DX42 

Difference between second 
table point greater than x 
and first table point less than 

X (X 1+ 2 - X i> 


46, 49 

DY21 

Difference between dependent table 
point corresponding to x. and the 

table point corresponding to 

vA - y ui ) 


57, 58 

DY32 

Difference between dependent 
table point corresponding to 

X i+1 an< * P°^ nt 


47, 48 



corresponding to x,, (y - y^ 



DYOX21 

Ratio of DY21 to DX21 
l(y i - y i.i ,/(x i - Vi" 


58-60 

DYOX32 

Ratio of DY32 to DX32, 

[ ( y i+i - y i )/( -Vi - x i>l 


48-50 

I 

Do-loop counter 


4, 5, 19-22, 86, 87 

1 11 

Subscript counter 


20-22 

IAR 

Dimension of the independent 
and dependent tables 


CALL, 2 ” " 

IFIRST 

First time entry indicator 


7, 31, 32 

IGO 

Internal loop control indicator 


28, 33, 35, 37, 
53, 68, 73, 76 
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TABLE 16. SUBROUTINE XNTERP NOMENCLATURE (Continued) 



I 

'I 


Symbol 

Description 

Units 

Reference 

IPOS 

Input array position indicator 


mmaarn 

IX 

Array subscript counter 


3, 9-13, 16-18, 
20, 25, 27, 30, 
40, 52, 63, 67, 
84 

IXIN 

First time entry indicator or 
saved value of subscript 
counter from previous entry 


CALL, 1, 84 

IX MAX 

Working maximum length of 
the input tables IAR - 1 


2, 8, 14, 18, 
30, 67 

IXO 

Internal limit for the subscript 
counter 


1, 6, 8, 25, 
40, 63 

IXOGO 

Subscript limit indicator 


26, 39, 51, 61, 
62, 85 

U1 

Ratio of DX2 to DX32 
t(x - x.)/(x. +1 - X.)] 


77-79, 81 

U2 

U1 squared 


78-81 

U3 

U1 cubed 


79, 80 

X 

Input value of the independent 
variable 


CALL, 11, 14, 
16, 23, 70 

XAR 

Input array of independent 
variables 


CALL, DIM, 
11, 14, 16, 21 

XI 

Array of table points surrounding 
the input value of the independent 
variable 


DIM, 21, 23, 
24, 46, 49, 55, 
56, 59, 70 

Y 

Output value of the dependent 
variable corresponding to the 
independent variable 


CALL, 82, 89, 
92 

YAR 

Input array of dependent 
variables 


CALL, DIM, 
14, 22 

YB1 

Value of dependent variable 
on the left parabola 


71, 82, 83, 89 

YB2 

Value of dependent variable 
on the right parabola 


74, 82, 83, 92 

YI 

Array of dependent table points 
corresponding to the XI values 


DIM, 22, 45, 
47, 49, 54, 57, 
59 


i 


208 





































TABLE 16. SUBROUTINE XNTERP NOMENCLATURE (Concluded) 


Symbol 

Description 

Units 

Reference 

YP 

Output value of the derivative 
of the dependent variable 


CALL, 83, 90, 
93 

YPB1 

Derivative of dependent variable 
from the left parabola 


72, 83, 90 


Derivative of dependent variable 
from the right parabola 


75, 83, 93 




SUBROUTINE ZETAIT 


Subroutine ZETAIT calculates the shape parameter £ = (6 /A) n 
and boundary layer thicknesses A, 6 , and 6 * at an axial point x for 
given values of 6 and <p . 

For £ s l 


where 



£ = 


(t h . 

\0 IJ + V£ 



For £ < 1 


6 * 

T 


( 


l/n - 1 R - I 7 
14 + I 5 , 


n(I 4 + I 5 ) ’ 




where 


1 



An iteration procedure is used to evaluate £ . 


(4) 

(5) 


( 6 ) 


(7) 

( 8 ) 


(9) 


1. An initial guess £ is made. 

g 

2. For the value of £ (- 1 or < l) appropriate functions of 

g 

P/P Q have been defined in subroutine FIIF, and ;he proper integrals 
are evaluated using subroutine INTZET. 



3. The term £ is calculated by the appropriate formula above. 


4. A relative error comparison is made. If 



TOLZET 


9 


convergence is considered to be satisfactory, If not converged, a 

form of Wegstein* s method (see step 42) is used to calculate a 

new guess for £ , and steps 2 to 4 are repeated up to a maximum 
S 

of 50 iterations. 


The boundary layer thicknesses 6*, 6, and A are then 
calculated according to equations (4), (5), (7), (8), and 
A =£ n 6. 


COMMON BLOCKS 


COMMON blocks COFIIF, INPUT, INTER, OUTPUT, and SAVED are 
used. 


TBL SUBROUTINES 

Subroutine BARPRO calls ZETAIT*. 
ZETAIT calls INTZET. 

FORTRAN SYSTEM ROUTINE 


No FORTRAN library routines are used. 

Built-in FORTRAN function ABS is used. 

CALLING SEQUENCE 

The subroutine calling sequence is 

CALL ZETAIT 
SOLUTION METHOD 

Save the value of <p/Q 

1. ERASE 1 = (PHI) /(THETA) 
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Set indicator to evaluate integrals... 

2. IFINT = 1 

An iteration procedure is executed to evaluate £ = (A/6 )^ /w i n 

the following Do-loop, using an initial guess £ . 

g 

3. Do 43 I = 1, 50 

Save n = (MZETA) used in the exponent of the velocity, equation 

4. MMINT = MZETA 

Save the wall enthalpy H = (A) defined in subroutine BARPRO. 

w 

5. AFINT = A 

Approximation of £ , initially using £ - (0/0) Ax. , calculated in 

g initial 

su^.outine BARCON. 

6. ZETAG = ZETA 

Check whether the temperature thickness A exceeds velocity 
thickness 6 . 

7. If ZETA ^ 1.0, go to 19 
If ZETA < 1.0, go to 8 

Save H - H = (B) which is determined in BARPRO. 
o w 

8. BFENT = B 

Calculate the dynamic enthalpy multiplied by square of £ with minus 
sign. 

9. CFINT = (C( (ZETA) 2 

10. CALL INTZET4 0.0, 1.0, ZI1P) 

Divide B = H - H by £ . 
o w 




21. CALL INTZET( 0 . 0 , 1.0, ZI1). 

Save H - H » (B). 
o w 

22. BFINT = B 

U 2 

(3 o 

Compute - -T- '(/, 
d* 

23. CFINT = (C)(ZETA) 2 
Set (i/e). 

24. ERASE2 = (1.0)/(ZETA) 

Ut a n 

Call subroutine INTZET to obtain l£ = / W (l - W)dW. 

o P e 

25. CALL INTZET( 0.0. EBASE2. ZI2P) . 

Save (n - l). 

26. MMINT « MZETAM 
Calculate (H - U 2 /2). 

27. AFINT = A + C 
Set 

28. CFINT = 0.0 

1 . 

Call subroutine INTZET to obtain I’ 3 = f — W-“-(l - W)dW. 

i/e P e 

29. CALL INTZET(ERASE2, 1.0, Z13P ) 



Calculate £ by an appropriate formula £ = 


n+1 


t 

B 



30. ZETA = ( ( ERASE!) ( ZIl) /( ( ZI2P + ZI3P)/(ZETA))) 

Determine the relative error (£ - £ )/£ . 

‘ S R 

31. DZETA = (ZETA - ZETAG)/(ZETAG) 

Check whether convergence is obtained. 

32. If | DZETA | < TOLZET, go to 46 
If | DZETA | 2= TOLZET, go to 33 

Check the Do-loop counter. 

33. If I ^ 2, go to 37 
If I < 2, go to 34 

Save £ in the case of I = 1. 

34. Z4 = ZETA 

Save £ . 
g 

35. Z2 = ZETAG 

36. Go to 43 


RMZETA 


Save Z4 and Z2 when I ^ 2. 


Set 


39. Z4 - ZETA 

40. Z2 = ZETAG 

41. Z5 = (Z4 - Z3)/(Z2 - Zl) 

Approximate a new £ using a form of Wegstein’ s method. 

g 

42. ZETA = (24 - (Z5) (Z2))/(i. 0 - Z5) 

43. Continue 

Print out the error indication of shape parameter iteration failure. 

44. WRITE X, ZME, THETA, PHI 

45. WRITE Zl, Z2, ZETA, Z3, Z4 

Proceed with calculation when convergence is obtained. 

Set IFINT equal to two. 

46. IFINT = 2 

Save n in the exponent of velocity profile. 

47. MMINX-= MZETA 

Save wall enthalpy H . 

w 

48. AFINT = A 

Save (H - H )/£. 
o w 

49. BFINT = (B)/(ZETA) 

Save ( -UV2) . 

50. CFINT = C 
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Calculate £ n . 

51. ZETATM = (ZETA) ZMZETA 
Check whether £ exceeds one. 

52. If ZETA 5: 1.0, go to 60 
If ZETA < 1.0, goto 53 

Call subroutine INTZET to obtain I 6 = J — S dS 
in the case £ < 1.0. o ^e 

53. CALL INTZET( 0.0, ZETA, ZI6) . 

Save stagnation enthalpy H^. 

54. AFINT = A + B 
Set 

55. BFINT =0.0 

1 

Call subroutine INTZET to obtain I 7 = f — S dS. 

£ P e 

56. CALL INTZET( ZETA, 1.0, ZI7) . 

Calculate (I 4 + I 5 ) . 

57. ERASE2 = ZI4 + ZI5 

I - Is - It 

Compute 6*/0 - T in the case of £ — 1, i.e. , 6 — A. 

a 4 + A S 

58. DELSOT = (OOMZET - ZI6 - ZI7) /(ERASE2) 

Q 

Obtain the velocity thickness 6 = — tt — ? , 

n (I 4 + I 5 ) 
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59. Go to 67 


1 

Call subroutine INTZET to obtain I 2 = / “ S n dS 
in the case of £ > 1, i, e. , £ < A. o e 

60. CALL INTZET( 0.0, 1.0, ZI2) . 

Save (n - l). 

61. MMINT = MZ ETA M 

Save (H - U 2 /2) . 
w e 

62. AFINT = A + C 
Set 

63. CFINT * 0.0 


Call subroutine INTZET to obtain I 3 = J S n_1 dS. 


1 P e 

64. CALL INTZET( 1.0, ZETA, ZI3) . 

B e 

Obtain the velocity thickness o = — r- 

nli 

65. DELTA = (THETA) /(ZMZETA) /(ZIl) 

, 5* £ n /n -1,-1, 

Compute — = ■ 2 — 1 " 

e ll 

66. DELSOT = ((ZETATM) /(ZMZETA) - Z13 - ZI2)/(ZIl) 
Obtain the temperature thickness A = £ n 6 . 

67. BDELTA = (ZETATM) (DELTA) 

Obtain the displacement thickness 6*. 
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68. DELSTR = (THETA) ( DELSOT) 


69. Return 


Subroutine ZETAIT nomenclature is given in Table 17. 


TABLE 17. SUBROUTINE ZETAIT NOMENCLATURE 


Symbol 

Description 

Units 

Reference 

A 

Wall enthalpy H defined 
w 

in subroutine BARPRO 

ft 2 /sec 2 

/SAVED/, 5, 14, 27, 
48, 54, 62 

AFINT 

Saved value of A 

ft 2 /sec 2 

/COFIIF/, 5, 14, 27, 
48, 54, 62 

B 

Stagnation enthalpy 
minus wall enthalpy 
calculated in BARPRO 

ft?/sec 2 

/SAVED/, 8, 11, 14, 
19, 22, 49, 54 

BDELTA 

Temperature thickness A 

ft 

/OUTPUT/, 67 

BFINT 

Saved value of B or 
B/ZETA 

ft 2 /sec 2 

/COFIIF/, 8, 11, 15, 
19, 22, 49, 55 

C 

Dynamic enthalpy with 
minus sign (-U^/2) 

calculated in BARPRO 

ft 2 /sec 2 

/SAVED/, 9, 12, 20, 
23, 27, 50, 62 

CFINT 

Saved value of zero 

(-U 2 /2) or (- £ 2 U 2 /2) 
6 6 

ft 2 /sec 2 

/COFIIF/, 9, 12, 20, 
23, 28, 50, 63 

DELSOT 

Displacement thickness 
divided by momentum 
thickness 8*/d 


/OUTPUT/, 58, 66, 
68 

DELSTR 

Displacement thickness 
6* 

ft 

/OUTPUT/, 68 

DELTA 

Velocity thickness 6 

ft 

/OUTPUT/, 65, 67 

DZETA 

Rfilative error 

(c - £ JA 

g g 

1 

31, 32 

ERASE 1 

Energy thickness divided 
by momentum thickness 

f/e 


1, 17, 30 

ERASE2 

Temporary storage 
variable 

— 

24, 25, 29, 57, 58 

I 

Do-loop counter 

— 

3, 33 
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TABLE 17. SUBROUTINE ZETAIT NOMENCLATURE (Continued) 


Symbol 

Description 

Units 

Reference 

IFINT 

Switch indicator for 
calculating integrands in 
subroutine FIIF 

— 

/COFIIF/, 2, 46 

MMINT 

Exponent for integral 
evaluation 

' 

/COFIIF/, 4, 26, 
47, 61 

MZETA 

n. in the exponent of the 
velocity profile 

— 

/INPUT/, 4, 47 

MZETAM 

n - 1 

— 

/INTER/, 26, 61 1 

OOMZET 

1/n 

— 

/INTER/, 58 

PHI 

Energy thickness 0 

ft 

/OUTPUT/, 1, 44 1 

RMZETA 

l/(n + 1) 

— 

/INTER/. 17, 30 

THETA 

Momentum thickness 9 

ft 

/OUTPUT/, 1, 44, 
65, 68 

TOLZET 

Tolerance for £ iteration 

— 

/INPUT/, 32 

X 

Axial distance 

ft 

/OUTPUT/, 44 

N 

Previous value of £ 

g 

— 

/OUTPUT/, 38, 41, 
45 

Z2 

Saved value of £ 

g 

— 

/OUTPUT/, 35, 38, 
40, 41, 42, 45 

£3 

Previous value of £ 

— 

/OUTPUT/, 37, 41, 
45 

Z4 

Saved value of £ 

— 

/OUTPUT/, 34, 37, 
39, 41, 42, 45 

Z5 


— 

/OUTPUT/, 41, 42 

ZETA 

Shape parameter 
£ = (8/A) l/n 


/OUTPUT/, 6, 7, 9, 
11, 13, 16, 17, 19, 
23, 24, 30, 31, 34, 
39, 42, 45, 49, 51- 
53, 56, 64 

2ETAG 

Approximation for £ 

— 

6, 31, 35, 40 

ZETATM 


— 

/INTER/, 51, 66, 67 

ZI1 

Il= f ~ s n (i - S)dS 
o M e 

— 

/SAVED/, 21, 30, 
65 t 66 
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TABLE 17. SUBROUTINE ZETAIT NOMENCLATURE (Continued) 


Symbol 

Description 

ZI2 

l = f JL 8 n dS 

1 J p 

O C 

ZI3 

1 3 = f — S n-1 dS 
o P e 

ZI4 

h = 1 -^-S n (l-S)dS 
o 

ZI5 

h = f ~~~ s n (i - S)dS 
t p e 


Units 



Refei’ence 


/SAVED/, 60, 66 



/SAVED/, 64, 66 



/SAVED/, 13, 17, 57 



/SAVED/, 16, 17, 57 



/SAVED/, 53, 58 


1 1 - ; f s* ds 

t e 


/SAVED/, 56, 58 


ZI1P 


I \ m f — W n (l - W)dW 
o P e 


/SAVED/, 10, 17 


ZI2P 


f — W n (l - W)dW 
o P e 


/SAVED/, 25, 30 


1 Vw 

1 I 
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TABLE 17. SUBROUTINE ZETAIT NOMENCLATURE (Concluded) 


Symbol 


ZI3P 


ZME 


ZMZETA 


Description 


1 5 - / — w n-1 (i - W)dW 

l/E P o 


Mach number in free 
stream 


Real value of n 


Units 


Reference 


/SAVED/, 29, 
30 


/OUTPUT/, 44 


/INTER/, 51 














APPENDIX. DERIVATION OF EQUATIONS [1] 


This appendix is a complement to Reference 1 and covers the deriva- 
tion of important equations in detail from section 2.2 to 2. 9. The equation 
numbers are identical to the ones in Reference 1. 

Derivation of Equation (7)(Displacement Thickness) 

The potential nozzle flow can be considered identical to the flow in a 

th 

real nozzle from the nozzle axis to the n streamline close to the nozzle 

th 

wall. In addition the flow field would extend a distance 6' from the n 

P 

streamline to a fictitious potential wall if the flow would possess the core 

flow properties. The real flow, however, extends a distance <5' from the 

th ^ 

n streamline because of the existence of a boundary layer. The mass flow 

rate rii through the potential flow-field extension 6’ must be equal to the 

P 

flow-field extension 6^ considering real-flow behavior. 


ih = m 
P r 


This postulation provides the key relating the real-flow condition to a potential- 
flow equivalent and makes it possible to determine the displacement 

thickness 6 which identifies the dislocation of the potential wall to account 
for the real-flow behavior. Introducing equations (1J and (4) into the previous 
equation results in 


27rr p U 6* 

00 00 p 


6 * 

/r 

pu dy 
o 


6 » 

P 


P U 

OO 00 



o 


pu dy 


» 


223 


6 ' 

r 


6 « 

r 


6' = 
P 


6' 

f r p u 
J p U 

O oo oo 


dy 


h 


6' 


■ / r * -J 


6 ' 

r 


_£iL_ 

p U 

oo oo 


■:y 


so that one obtains the displacement thickness 6* 



n* STREAMLINE 



Derivation of Equation (10) (Momentum Thickness) 

Using equations (2) and (5) , the deficiency of momentum flux in the 
real flow as compared to the potential flow can be represented by 
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Si - Si = 2irr p U 2 6' - 27rr f r pu 2 dy 

P r 'oQ OO p J 


= 2nr p U 2 6' 


(«• - / r dy) 

V p o p U 2 / 


Introducing equation (71 yields 


6- = 6' - f r ( 1 - -CSL W , 

P r ”o V P co U oo/ 


Si - Si = 27rr p U 2 d' 
p r oo oo r 


d» / v 6* „ 

-rf'-fi-)*-' 

o \ oo oo / op 


= 2?rr p U 2 / r - 

oo oo J n U 
O » oo 


(■•t) 


dy 


where the integral represents the momentum thickness 6 . 

Derivation of Equation (13) (Energy Thickness) 

Applying equations (3) and (6) yields the deficiency of enthalpy flux of 
the real flow as compared to the potential flow. 
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I • 


6 ' 


E - E « 2nr p U (H - H ) 6» - 2nr f p u(h - H ) dy 
P r « oo o w p J H ' o w 7 y 

o 


27rr p U (H - H ) 

00 00 o 


6' / h - H 

6’ - f r 

p J p u \h - 

o oo CJD \ o 


W 


H 


dy 


w 


Using equation (7) results in 


6 ' 


% - 5 r - / r ‘ - fr * ' 

O \ oo OG / 


E - E = 2fr p U (H - H ) 
p r w oo o w 


p 

6' 

/ v 

6 ! - 

I r 

i . z ) 

r 

J 

\ P u 


o 

\ oo oo y 


dy 


6' 

h 

- H \ 

* 

r r pu / 

o 

w ' 

1 dy 

P u 1 

y H 

- H / 

0 oo oo 

\ o 

w / 



= 27rr p U (H - H ) 
oo oo o W 


f 6 r pu / h o ~ H W \ 

J pu l ~H - H ) dy 

O oo oo \ 0 W / 


The integral represents the energy thickness cp , 

Definition of Boundary Layer Thicknesses 


Displacement Thickness 
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Momentum Thickness 


Q = 


/ 


6 ' 


r _£U_ 

P U ' 


Q oo (O 


u 

U 


idy 


Energy Thickness 


O oo oo \ 


- u dy • 


Derivation of Equation (16) (Skin Friction Coefficient) 

The friction coefficient is defined by the force in direction along the 
F 

wall per unit area — = r made dimensionless by the dynamic pressure 
Aw 


2r 

C = — — . 

1 P u 2 

OO 00 


Equation (17) (Stanton Number) 

This dimensionless parameter can be interpreted as the ratio of the 
heat flow density emerging from a wall for each degree of temperature 
difference between the wall and the fluid and the heat flow density convected 
by the flowing medium for each degree of temperature difference. 


St 


Nu 

EePr * 
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hg = specific heat transfer coefficient 

(contains all influences caused by properties and flow conditions) 

\ = heat conduction coefficient ( is a material constant) , 

I = characteristic length 

Therefore, the Stanton number can be expressed by 


St 



> 


JL 


C p u 

p «0 00 


I 


C = specific heat at constant pressure, 



The total heat flow between a surface and a fluid is 


Q = h A (T - T ) 
g F w 


•f- - q -- h (T - T ) 
A w g F w 


h g = 

F w 


Substituting this expression in the Stanton number yields 


C p U (T 

p oo oo jr 


T ) 
w 


or since 


h = C T for C = constant, 
P P 


p U (H - H ) 
oo oo aw w 


Since T is the temperature which a specific location on the surface assumes 
F 

when the convective heat transfer is zero, it can be considered as the recovery 
temperature. The corresponding enthalpy expression is the adiabatic wall 
enthalpy. 


Equation ( 18 ) (Adiabatic Wall Enthalpy) 


When a gas flows past an insulated or adiabatic surface, the tempera- 
ture at the surface will rise above the temperature of the free- stream gas 
but will not quite reach the total temperature. The temperature at an adia- 
batic surface is called the adiabatic wall temperature. 

In practice it has been found convenient to relate T and T by 

aw °o 

the recovery factor R t , which is a measure of the fraction of the free- 

stream, dynamic-temperature rise recovered at the wall. The term R 
is defined as ” 



T - T 

aw °° 


U 2 



Jg 


2 J 

(y - 11 M 2 

' OC 



The derivation of this equation is shown below. 

T aw - T » = * T aw - 2 ° D * J 

U* , a ! 

00 U 

2c gJ 00 a 2 

P 

U 

where speed of sound is a = n/v g R T ' Mach number M = — 

* a a 

00 

Therefore, 


T - T 

aw oo 


U 2 


(T 


aw 


M 2 

a 


- T ) 2c gJ 

P 

ygRT^ 


Introducing the thermodynamic relationship 



C 

v 


results in 


T - T 
aw 00 


U 2 


2 c g J 
P 


T - T 
aw 

T 


2 c J 
P 


y (C - C ) M 2 
r p v a 

* r* 




OO 



TEMPERATURE PROFILE 
FOR HEAT FLOW FROM 
THE SURFACE 





For turbulent boundary layers the recovery factor R can be approximated 
by the relation ^ 



3 

~ ^ p r ' 


The ratio of the adiabatic wall enthalpy and the stagnation enthalpy is defined 
as 



h + R 


U 2 
00 

T T~ 


u 2 

h+ — 
2 


i 


where R^ — recovery factor [ the fraction of dynamic enthalpy ( or tempera- 
ture) which, is recovered] , 


U 2 

OO 

~ 2 ~ ~ ^ T 0 " V ~ dynamic enthalpy, 

h = free -stream static enthalpy. 

Rearranging the previous equation yields 

U 2 

OO 

T 2h 
U 2 

OO 

T h 

and introduv ing the recovery factor R = Pr ^ for the turbulent boundary 
layer results in 
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H 

aw 

H 


1 + Pr 




2h 


1 + 


U 2 

oo 

~2h 


Derivation of Equation (21) 


3 




Since the gradient of momentum flux for the potential stream tube of 

area A = 27rr 6' should be balanced by the pressure gradient acting over 
P 

the same area A , we obtain the following equation. 


/ 

( ft + — r 5 ax 
V P dx 


ft = p - ( p + 4^ ^ x ) 27rr • 

p L dx /J p 


Thus 



dP 

- 27rr 6' ~r 
p dx 



3 From now on the subscript °° is omitted for convenience. 


I 


Derivation of Equation (23) 

For llio real flow the viscosity effect appears, 
equation in the previous section is modified as 

/ dM x / 

(V d/ AK ) V [p- Ax) 

Thus 


dM 


r .. dP 

f .~ - - ^ /ri ’ ^ , — - 27rr r 
dx r dx w 


Derivation of Equation (24) 

Beginning with equations (22) and (23) 

d „ * d P 

77 (M + 27rr pU 2 6 ) = - 2irr (6' - 6 ) 7^ 

CIX L r d x 

(iCl > - 27rrT - 27rr6' 

dx r w r dx 

and subtracting these equations from each other yield 

~ (27rr pU 2 6 ) = 27rrT + 27rrd" 7 ^ 

ax w dx 

After factoring out 2 tt , one obtains 


do h p 

T~ (rpU-0) =■- r r + r<5 — - 
dx w dx 


Therefore the first 


2?rr 6' - 27rr Ax r 
r w 
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6 


d£ 

dx 


£ 

dx 



e 


/ 6 * 

(£_ilH + i^ + IdU + ldp 1 
\ U dx U dx U dx p dx "* r 

/ 6 * 

( - - 1 + 

\ U dx U dx p dx r dx 


Consider only the two terms 


^1 dU + 1_ dp 
U dx p dx 


According to the mathematical definition 


d In x _ _1 
dx x 

one can write 

_1 dU 1 
U dx + p 


dp _ d In U dU d In p _dp 
dx dU dx dp dx 

d In U + d In p 
dx dx 

- x ( 1" 0 + In P) 


• to 1 ln (urt) 

_ d [in ( Up)] dUp 
dUp dx 

_ _i_ dpu 
“ pu dx • 
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Substituting this into the previous equation results in 


C / — + 1 

d fi 2L i du 1 dpU | 1 dr 

dx 2 ^ V U dx ' pi) dx ' r dx 


Derivation of Equations (26) and (27) 


Geo me try 



Substituting this into equation (251 yields 



This equation is applicable for axisymmetric flow. 
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For two-dimensional flow 


^ and 


l 

v 


dy 


dx 



dz 


dx =:: [(dy) 2 + (dz) 2 ]^ 

% 


dx 

dx. 


d/. 


+ 1 


Substituting this into equation (2d) yields 


dO 

dz 


1 + 


(fj 


W 




- 0 


i 

o 


u 


1 + (f) 


dU 


72 dz 


pu 

“l + (^) r 



\ dz/ 



dpU 


d0 

dz 


‘ * (I) 


T 1 


h 


1 + — 

fl_ dU J_ dpU 

U dz + pU dz 


Derivation of Equation ( 30 ) 

The change in energy along a streamline for real flow is equal to the heat 
transferred between the fluid and the wall 

dE = 2?rrq dx 
r w 


dE 

i 

dx 


= 27rrq 


w 


( The sign of the right-hand term depends on the direction 


ion of q . ) 
w 


Heat transfer from the wail to the fluid is indicated by a negative sign 
Another way of the derivation is given below. 
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/ clE \ 
E r -{K * “* J = 


dE 


= 27rrc 


dx - r \r 

Derivation of Equations (32) through (34) 

Beginning with equation (3l) and differentiating by parts yield 


(30) 


y- [rpU (H - H ) <p] = rq, 
dx o w 


w 


rpU (H - H w ) g + rpU* 


<i(H - H) 


— + rp(H - H )<£ ” 

dx o w dx 


dp dr 

♦ rU(H o • V* to * PU(H 0 ' V* to = rq w 


Substituting q w by equation (17) 


q = C„pU(H - H ) 
H aw w 


ar.d solving for ^ yield 


Ww - V - ^ 

dx 


d < H o - v 

(lx 


- rp(H - H )4> ^ - rU(H o - HJO ^ - pU(H Q - HJ<t> ^ 


r P U(H o - V 
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. . h - h 

d0 _ aw w 

dx C H H - H ^ 


o 


w 


H - H 

o w 


d(H - H ) , , 

£ w_ 1_ dU 1 dp 

dx U dx p dx, 


1 dr 
r dx 


d(H 


H ) 
w 


dx 


dll 
c 

dx 


dH 


w 


dx 


Since H is assumed to be constant, 
o ’ 


dH 


dx ° ; 

1_ dU 1 dp _ 1 dpU 

U dx + p dx ' pU dT ' L See derivation of equation (25).] 

Substituting these terms into the previous equation yields 

1 


j . H - H 

d£ _ aw w 

dx C H H - H 


<P 


w 


H - H 

o w 


dH « . tt 

w 1 dpU 

dx + pU dx + 


1 dr \ 
r dx J ' 


(32) 


For axisymmetric flow 
dx = (dr 2 + dz 2 )^ 

Vt 


** - [(£) ! + ‘] i/! dz ■ 


This equation substituted into the previous equation yields 


ha / H - H 

dg _ „ J aw w 

dz H \ H - H 
o w 


1 + 


V 2 


- <p 


H - H 

o w 


dH 

w 

dx 


+ M + i £ | m) 

pU dz r dz ' * (33 * 


m 
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For two-dimensional planar flow r 


CO 



This equation substituted into equation (32) yields 




J_ dpU \ 
pU dz j ' 


(34) 


Derivation of Equations (36) through (39) 


The skin friction coefficient in a nozzle is assumed to be the same as 
for a flat plate at the same free-stream conditions p , U , p , T , M , the same 
wall temperature T w , and the same momentum thickness 6 . Unfortunately, 

even this drastic assumption does not permit a completely reliable evaluation 

of C„ since only the adiabatic skin friction coefficient C. , obtained when 
f fa ’ 

T = T , is known accurateiy. The relationship between C. and C„ for 
w wa r f fa 

severely cooled turbulent boundary layers, when gas properties vary greatly 

between the free stream and the wall, is uncertain. 

The most accurate method for computing the adiabatic skin friction 
coefficient is judged to be that of Coles [ 11] . He employs a transformation 

method which gives C as a function of Mach number M and the Reynolds 

l2L 

number based upon the momentum thickness as defined by 
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For largo temperature differences the properties of gases cannot be 
assumed constant. They depend very much on the temperature. Using the 
Prandtl number, which contains basically only the property parameters of a 
fluid, 


Pr 


v_ 

a' 


a' 


A 

V 


Pr 


vC p pC 
P P 


A A 


Since the dynamic viscosity p and the conductivity A vary considerably 
faster with temperature, one can rearrange 

p _ Pr 
A “ C 

P 

This equation shows that for a constant Prandtl number Pr and constant 
specific heat C , the heat conductivity A varies proportionally to the viscosity 

p . Therefore the law of temperature dependence has to be prescribed for one 
property only 

P ~ T 

For ideal gases Reference 12 gives an equation for the dynamic viscosity 

p = Kn/STCT 1 K is a constant 

91T = molecular weight 


or 


Assuming a reference condition (01 the above equations yield 



P T 
o o 


( 36 ) 


The exponent m is dependent on the temperature and varies between 1, 0 
and 0. 5. 




It is interesting to note that for m = 1 the dependence of the friction 
factor on the Mach number and the ratio of T (wall temperature) and 

(stream temperature) vanishes. It may therefore be expected that property 

values introduced at a proper reference temperature situated somewhere 

between the temperature extremes, encountered within the boundary layer, will 

cause the variation of the friction factor with M and T /T to vanish. It 

aw *> 

has been shown that such a reference temperature exists. Using the latter 
relationship 



the equation for 


can be rewritten 



In this equation the term 
as derived below. 



and 


(pU) 


can be replaced by relationships 


Using the equation derived in step 7 yields 



<30 


The expression pU is derived below. Starting with the speed of sound relation- 
ship 




and the ideal gas equation 








■ t 




f - RT 

results in 

a 2 ~ 2JL . 

P 

Applying the Mach number definition 

M = — or M 2 = , 

a a 1 ’ 


M 2 = 

yP 


(multiply and divide by p ) 


1 


M 2 = 


ypp 


P 2 U 2 = yppM z 


>k 


pU = M (ypp / 2 

Modification of the term in parentheses results in 


ypp = ypp 


P P 
o o 


P P 
oo 


pp 

= y — *— p p , 
1 p p o o 
oo 


— = RT 
P o 


P o RT 
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i! 


V 

I 





together with the adiabatic relation of temperature and pressure 


7-1 



and the previously derived equation 


T 1 + 2 


r-z-1 M 2 


yields 


7-1 


(?) 


1 + y ~ 1 M 2 

dt 


7 ■ (‘ + ^ “*) 


7 

7-1 


Substituting these terms yields 


7PP = 


7 2 P 2 



y+1 

7 


yRT V p 
' o ' o> 


? 2p o 2 

■ ‘ilo. 


ll 


(‘ - m! ) y 


y+1 

7 


Introducing this expression in the original equation (page 244) yields 



Coles (11| shows that the available data on the variation of C 

fa 

(adiabatic skin friction coefficient) with R^ and M_ can be correlated within a 
few percent by a single curve 



_V_ C 

M p fa 
aw aw 


Multiplying by R^ yields 


(37) 


C A ■ 


PM 

P M 
aw aw 


C, 
ia e 


(38) 


where p and u are the density and viscosity evaluated at T and u 
aw aw aw s 

is the viscosity at a mean sublayer temperature T given by 

s 


T 

s 


T 

aw 


+ 17.2 





. (39) 


It is evident that C f and R are merely values of C, and R for low-speed 
i 0 I a. u 

flow. Coles' relation between C, and C„R is shown in the Figure A-l. 

lie 



Figure A-l. Coles' relation between C^. and C^R^ . 
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For values ol' C^lt^ above 64.8, the relation employed is the asympotic one 
given by Coles l llj . 


2. 44 In 


c r, 781 . 

L f ' * ' 1 / 2 y / 2 




+ 7.68 


To permit the choice of a small initial boundary layer thickness the 
range below C^R^ = 2.51 is extended by 

0.0098 96 

C f = (CR J 0 - 562 
1 6 

Writing density and viscosity in terms of temperature, using the ideal 
gas equation 

p 

— = RT 
P 


and considering that at the same location P = P , one obtains 

aw 


T 


Applying the previously derived relationship of viscosity and temperature yields 


M P 

M P 
aw aw 


_ \m T 
T \ aw 


/ T aw 

V T 




Introducing the proper terms in equations (.37) and (.38) results in 


'fa 


(37a) 


(V)(W“ 

aw 


/T \ Um 

V*. - (-f) 


C ta ll „ 


(38a) 


Thus given all the gas properties outside the boundary layer and given 6 , 
then can be calculated from equation (35) . A value of is then assumed 

and equation (38) is used to find C^R^ . However, C^R^ is related to 
by the data in Figure A-l and equations (37a) and (38a), Therefore, can 
be found from C^R ; but is related to by equation (37). Therefore, 
the calculated is compared to the assumed and an iteration is per- 
formed until convergence is obtained. The skin friction coefficient is then set 
equal to , which means that there is no effect of heat transfer on the skin 

friction coefficient. 


Derivation of Equation (40) 

It is assumed that the Stanton number C u in a nozzle is the same as 

on the flat plate at the same free-stream conditions p , U , j it , T , M ; the 
same wall temperature T w ; and the same energy and momentum thickness <p 

and 6 . The most accurate available relation for computing the flat plate 
Stanton number is von Kar man's form of Reynolds' analogy which relates 
Stanton number to skin friction coefficient, with a secondary correction for 
non unity Prandtl number as follows [13] . 
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This equation is valid when the momentum thickness and energy thickness are 

equal ,? - 0 . Under the same circumstances R =11, , where R , is the 

0 cp 0 

Reynolds number based upon the energy thickness defined by 


R 

0 fi 


If R^ is employed instead of It^ in the procedure for computing the skin 

friction coefficient, a number is obtained which is designated Cj(ll ). If 

0=0, then CRRj ) = C, and the Stanton number in this special case is 
10 1 


C 


H 


C f (R ) 
£ £ 
2 


r c f<y] 

1 r 

2 




Since C must depend more on <p than on 0, this equation is considered 
H 

as the first approximation to for unequal 0 and 0 as well. (The term 
C„ must approach infinity as 0 approaches zero, regardless of the value of 

0 . ) When test data were compared with the previous equation, it was found 

that the Stanton number and skin friction coefficient are insensitive to T /T 

aw o 

for cooling, which was assumed when . 


The effect of nonunity 0/0 cannot be large for 0< 0 . However, test 
data show same scatter, and a larger effect could be present for nozzles, 
since here <p >0, To include such an effect, the final equation for C is 

n H 

multiplied by (0/0) allowing for a correction. 



( 40 ) 



While C (R ) varies approximately with the -*/ power of 9 (Blasius 

f 9 

equation), the last equation for contains the following approximate 
dependence on 9 and 0. 1 


C 

H 


1 / o> ~ 

U -n .n 

e e 


Since C must depend more on 9 than on 0 and, in particular, must 
approach infinity as 0 — 0, it is evident that the upper limit for n is n 
An estimate of the actual value of n was made with a relationship of 


= % 



which is the same as the previous equation for n - 3/28 ^ 0. 1 . Thus n 
lies between zero and 0. 25 with some basis to choose 0.1. 


Derivation of Equations (43) and (44) 

Q = q A = 2JTr x q , 
w 4 w w 

dQ = 27rr q dx 
w w 

For axisymmetric flow 
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(43) 


«,[(§)’ + *] ! dz • 
For two-dimensional planar flow 



Therefore dividing by 2?rr yields the heat transfer per unit width 



Derivation of Equations (45) through (47) 

The total drag along a surface is represented by 


(44) 


D ■ Cfi £• . 

For axisymmetric flow the drag in z-direction can be determined if the 
. angle between the wall and the z-direction is known. 
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1 + tair a 


cos a 



A = 2m' x , 

D = 2 r7rC f ^ U 2 x , 

The force in z-dircction is 

F = F cos a , 
z x 

dD = 2*r C f | U 2 cos o dx 



dD = 27rr C f | U 2 

" / , \ 2 
(df) + 1 

*4 

cos a d 

z 

D = / 27rr C f - U 2 
o 

/ drV + j 
\dz / 

v 2 

COS Q! dz 


For two-dimensional planar flow 

I* ^ 

27rr 00 


Dividing the equation by 2?rr the drag per unit width is determined 
Furthermore 












Geometrical Representation of Possible Velocity and 
Temperature Thickness Combinations for the Boundary 
Layer and the Derivation of Appropriate Equations 



ij; . 

i ' 


g(y> 




I. 6 < A 

(a) y < A 

(b) 6 < y < A 



(a) y < A 

(b) A < y ^ <5 



Case I 

6 ^ A : The velocity boundary layer is smaller than or equal to the 
temperature (enthalpy) boundary layer. 

(a) y < A In this case y can be within or outside of the velocity 
boundary layer. 


i 



’ | 

. I s 
£ 

* 


I 
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For the whole range the velocity boundary layer is 

/ \l/n 

&-(*) 


Ify>6-»S = 1 Outside velocity boundary layer. 
Ify<6-*-S<l Inside velocity boundary layer. 

For the whole range the enthalpy (or temperature) boundary layer is 


h - H 

o w _ 

H - H 

0 w 


■(«’ 


= W ^ 1 


Ify<6-*W<1 Inside velocity boundary layer. 

If S < y — A-*-w^ 1 Outside velocity boundary layer. 

When the velocity and enthalpy ( temperature) boundary layers are 
compared, 


M 

(D 


1 7n W ^ 


and the following definitions result. 
1/n 

s = (f) ■ 


-a) 

■ (ff 


Therefore the velocity boundary layer can be expressed by 
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*1 NH #>- 


u = U S for S < 1 ; 

whereas the enthalpy (temperature) boundary layer yields 


(59) 


h - H 
o w_ 

H - H 
0 w 



(60) 


h = 
o 





( 61 ) 


with h Q the stagnation enthalpy within the velocity boundary layer. The 
stagnation enthalpy is composed of a static and dynamic portion. 


h = h + £ 2 . (62) 

O 2 

Substituting this in the last equation and making use of u = U S equation 
(65) is obtained. 


h = H + f (H - H) - . (65) 

w f, o w 2 

The density variation can now be determined from the universal gas 
equation 

— = Rt (Pressure is constant across the boundary layer.) 

P 



where t, temperature in the boundary layer, can be calculated from h- by 
means of 

h = C t 
P 
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However the results are only valid for 


W 


(b) 6 < y < a In this case y is greater than the velocity boundary 
layer but smaller or equal to the temperature bound- 
ary layer. 

For this range the velocity boundary layer is 


l/n 

Ml) -■-* 


(69) 


The second integral of equation (75) is zero since it considers a portion 
outside the boundary layer, whereas the first integral is modified according to 


■(f) 


l/n 


S6 1//n = y 1//n 


_n. 

S 6 = y 


dy 


n-1, 


02- = nS - 6 
ds 


dy = nS 11_1 6 dS , 


and finally reduces to 

1 - 

e = n f £ S n (1 _ S) 6 ds 

*2 O 


(76) 


In order to obtain the equation for the displacement thickness 6*, the 
relationship used in the derivation of equation (76) is used. Furthermore the 
first integral on the right-hand side of equation (80) is evaluated. 



n f S n_1 ds 


n — - 0 

n 


and the result appears in equation (81) . In the third integral on the right-hand 
side of equation (80) the relationship 


± = 1 
U 


is applied since y is greater than the velocity boundary layer. 

To obtain the energy thickness in equation (85) , use is made of the 
relationship as de, eloped for equation (76) . Furthermore the following rela- 
tion is applied. 

S = w; , 

ds _ r 
dw " ^ 


ds = t, dw 


Case II 


6 ^ A 


(a) y 


The velocity boundary layer is equal to or greater 
than the temperature boundary layer. 

In this case y is smaller than or at the most equal 
to the temperature boundary layer. 


(b) A < y < 6 Here y is greater than the temperature boundary 
layer but smaller than or equal to the velocity 
boundary layer. 

Equation (10l) is equal to one because y is outside the temperature 
boundary layer. 




For the momentum equation (106), the relationship shown in connection 
with equation (76) is utilized. 

Also equation (ill) representing the displacement thickness 6* uses 
the relationship applied in equation (76) . The first integral on the right-hand 
side is solved and results in 

/ ns ”" 1 ds = S n = 1 . 

o o 

For the momentum equation leading finally to equation ( 116) , the 
second integral is equal to zero since it represents a quantity which is outside 
the temperature boundary layer. The first integral on the right-hand side is 
reduced by substitution of the relationships developed in equations (76) and 
( 85 ) resulting finally in equation ( 116 ) . 


Thrust Degradation 


In a real flow the boundary layer and its viscous effects cause an I 

SP 

degradation. Furthermore, the potential flow field is narrowed because of the 
presence of the boundary layer. 



Subsequently the derivation of the equation representing the degradation 
of thrust caused by the viscous flow effects in the boundary layer is presented 
in the following diagram, where p represents a contour of the potential (invis- 
cid) flow and r represents a contour of an exactly equivalent viscous flow such 
that the pressure P along the normal to the wall y is the same on r as on 
p and the mass flow through r - n is the same as through p - n . The 
velocity distribution is confined between n - r . 
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r (REAL- WALL) 


NOZZLE 

CONTOUR 



p (POTENTIAL 
WALL) 


_____ n th STREA MLINE 


a 


c • 8 ' cos a 

6p cos a ["""\ P 


Vr 


u cos a 


The mass flow rate of the potential flow through the area between p 
and n is 

i 

m - 27rr<5 Up . 

P P 

The mass flow rate of the viscous flow through an area between r and n is 

m = 27rr6’ Up 
r r 

However U and p are varying with y such that 
5 ' 

t /* £ 

6 Up = pu dy 

v» J 


and 


m = 27rr f r pu dy 
r o 
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The momentum flux in % -direction is 
fa = AiU 

\ ■ 

such that for the potential and viscous flow one obtains 
(fop) * = (2*r 6^ UpJ (U cos a ) 

O' 

(fo ) = 27rr J 1 pu 2 cos a dy 

o 

Because of the definition, a is assumed to be constant between n and r, 

6 ' 

r 

because — « 1; therefore, 

r * 


O' 

(fa r ) = 27rr cos a f r pu 2 dy 

o 

Since the same mass flows through n - p and r - n, one obtains 

ih = m , 

P r 

O' 

27rr0 ’ Up m = 2rrr / r pu dy 
P o 



Subtracting 0 ' on either side 



I I 

6 - 6 
r p 


S’ 

6’ 

r 

/» 1’ 

dy - / 

o 

O 

S' 


r 

o 

U - ^ ri 

V p u 


p U 


\ dy 


dy 


As defined earlier 


r 1 

6 "' =6 - 6 
r p 


’’ ■ f (■ - d) " 


Now, both momentum equations will be subtracted from each other 


/Ki > — (iCl i = 27rr<5 Up U cos a - 2nr cos a [ r pu 2 dy 

v p z v r z p o 


. 6 ’, "2 \ 
= 2nrp U 2 cos a ^ - f 1 dy J 


= 2nrp tr cos 


= 27rrp U 2 cos a 


■ «(v 6 " - f r jlr dy ) 


S’ 6’ . ... 

r * - r (» - fo) dy 

n O ' f 


6 ' 

- r fw dy 

o 


= 27rrp U 2 cos oi J 


5* 

r 


-(■• Sr )' 


pu 2 

ru * 


dy 


203 


6' 

= 2iTrp U 2 cos a J 1 
o 

6 * 

= 27rrp U 2 cos a J r 

o 

The integral expression however represents just the momentum thickness as 
derived in equation (12) . 

•- f v f\ h ( l ■ ! u) dy ■ 

o ' ' 

Therefore 

(l5l ) - (Id ) = 27rrp U 2 cos a 6 

' pz r z *■ 

Since the momentum flux is equal to the forces 

& = F 

and the forces F can be expressed by pressure P times area A 
F = PA 
one obtains 



151 = PA 




Differentiation yields for the z-direction 


dM _ A dP dA 

dz dz dz 


Assuming that the variation of the cross-sectional area at a station z to be 
very small and negligible results in 


d lil . dP 

dz dz 


f 

k 




264 


t 


For the potential flow 


A = 2rrr6 , 

P 


i * 

A = 27rr<5 = 27rr6 cos a 

z pz p 


For the viscous flow 


A = 27rr6 r , 


i » 

A = 27rr<5 = 27rr<5 cos a 

z rz r 

Furthermore the shear forces must be considered for the viscous flow with 

t 

t A - mu = $1 
w 

In this case the surface area along x must be considered 

t 

A = 27rrx . 

Differentiating the momentum equation yields 

dl& „ dx 

- — = t 2irr ~r~ . 

dx w dx 

dl5l 0 

- — = r 27rr 

dx w 

Using the previously derived relationships results in the following 
momentum equation in z -direction 

d • dP 

= - 27rr<5 cos a — 

dz p dz 

(A negative sign was chosen since the pressure decreases with increasing z.) 


d <«y z . dp 

=-=■ - - 27 rr <5 cos a — - 2irr t 

dz r dz w 
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Using- the previously derived relationships yields 

(iCl > -(Ml ~ 27rrp U 2 cos a 6 
p x ' r z 


i t .. 

6 - 6 m 6 ' 

r p 


for the potential ease results in 


(I (til ) 

p z _ d 


- [(til ) + 2irvp U 2 e cos al = _ 2irr (6 - 6') CO s a 

z l v z <*> J ' r 


Subtracting the equation representing the viscous case 


d (M ) 
r z = 

dz 

from this one yields 


-27rr6 T cos a. — - 2 nr r 
r dz w 


— * (27rrp U 2 e cos a) = 27rr6" cos Of ~ + 2*rr r 
ciz dz w 


Multiplying the equation by — results in 

27T 


d ( r p U 2 o cos a) = r<5 cos cv dP + r r dz . 

w 

Solving for the r term and integrating result in 

J* r T \v = ^ ( r P U 2 0 cos a) - J* r 6* cos a dP . 

According to a previous derivation, the force acting on the exit plane is for- 
1 . Potential flow 


(F ) = m U + P A + f ep P dA 

D Ze H n J Y 




2. Viscous flow with the inclusion of shear forces 



ze 


mU+PA + f er PdA - / er T dA’ 

n n J V * MJ T 


C C 


A 


w 


(The subscript c refers to chamber. ) 

The degradation in force caused by the viscosity and the shear forces 
is equal to the difference between the potential and the viscous flow 


(AF) = (F ) - (F ) 

v ' ze ' p ze r ze 


A 


I 


ep 


P dA 



+ 



r dA' 
W T 


The change in area normal to the z -direction as a function of z is 


dA = 27rr dr 


dA 

dz 


27rr 


dr 

dz 


dA = 2?rr ^ dz 
dz 



For potential flow 



For viscous flow 


/ \ d(r + 6* cos a ) 

dA = 2Tr(r + 6* cos ot) dz dz • 

For shear forces 

dA’ = 27r(r + 6* cos ot )dz 
T 

Introducing these expressions in the previous equation yields 

re dr re„ , ** \ d(r + 6* cos a) 

(AF) = / e P2Jrr — - dz - J P2ir(r + 5* cos ot ) ^ 


e o 


dz 


t-i 

+ f e r 2 tt ( r + 6* cos ot )dz 
J w 


= f e P27rr dz - 
J dz 


r z 


f e P27rr ~ dz + f e P27r6* cos a ~ dz 
J dz J dz 


Z e„o_ d(5* cos a) dz 


+ J e P'27rr 
o 


dz 


re _ d(6* cos a ) ^ 

+ J P2tt6* cos a ^ dz 


+ J e r 27r(r + 6* cos ot )dz . 
o 

The first two integrals cancel each other. 

. ra J\. dr d(6* cos a ) 

(AF) z = -J® 2ttP 6* cosa^+r ^ 

e o L 

CJ . d(6* cos Qt)L„ 

+ 6* cos a dz — dz 


tUk 

+ f e t 27r(r + 6* cos ot )dz . 

j w 
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In the last equation the following terms can be expressed differently. 


6* cos 


dr d(6* cos at) d , \ 

a — + r — ? = -7” (ro* cos a ) 

d z dz dz 

c 1(6* cos a) 1 f d o 2 \1 

5* c03 a -S— ^ 5 < 6t! “)J • 


such that 


c o 


~ (r6* cos a) + - t- (6* 2 cos 2 a ) 
dz 2 dz 


(AF) ^ = -J C 2irP 

+ f e 27 r r (r + 6* cos a)dz 
J w 


dz 


Since terms of - — C ° - S — « 1, they can be neglected especially when higher 


order terms are used. 


z z 

(AF) = -J e 2irp d (r 6* cos a ) + J e 2 T 7rr dz 
e o o 

The previously derived expressions on page 266 

/ r t w dz = J d (rp w U 2 0 cos a) - / r6* cos a dP — 

is substituted in the equation above, such that 

z z 

(AF) = -J 6 27 rp d (r6* cos o') + J 6 2/rd (rp ffl U 2 0 cos a ) 
Z e o o 

z 

-J e 27rr6* cos ce dP . 



Combining the following integrals yields 


f Pd (r6* cos a ) + f cos oi. <IP = f d (Pr6* cos o') 

The equation for (Af) can be written 

z 

e 

z z 

(AF) ^ - -2tt f c d (Pr6* cos rv ) -i 2ir f c d (rp^ 11*0 cos a) 
co o 

z z 

= -27rPr<5* cos a | + 2rrrp U 2 0 cos a I c . 

O o 

Since the displacement thickness 6* and the momentum thickness 6 
are zero at z = 0, the thrust degradation is obtained as 

(AF) = (_27rp r5* cos a + 2;rrp U 2 0 cos a) 

z e 0 exit 

or 


(AF) z = T 2nrp U 2 0 cos a ll - 

e \ / J exit 

The result is only a function of exit plane conditions. 

Derivation of Equation (122) 

First the term — — of equation (121) will be derived. The 

energy equation states that 




Multiplying both sides by y yields 



0 = dh + UdU 



Substituting the previously derived term yields 


0 


yll 


1 


dT + UdU 


Dividing the equation by a 2 ^ yHT yields 
0 


1 dT U 2 dU 
y - 1 T + aT U 


0 , _1_ dl’ + m 2 - (| y 

v - 1 T U 


with M 


£ 

a 


dT 

T 


(1 - yJM 


2 £!U 
U 


dU 


Now a relationship for — will be derived. Starting with the Mach 
number definition U 


M = — or M 2 = 


and representing the speed of sound a by 


a = \[yWT 
result in 

m«=- u! 


yRT ' 

Differentiating this equation in logarithmic form yields 

2 In M = 2 In U - In y - In R - In T , 

2 d( In M) = 2 d( In U) _ d(ln y) d(ln R) d( In T) 
dz dz dz dz ~ dz 
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d(ln M) dM _ d( In U) dU _ d(ln y) dy d( In R ) dR 
dM dz dU d z dy dz dR dz 

d( In T ) dT 
dT dz ’ 


2 dM _ 2 dU l dy 1 dR _ 1 dT 

M dz ~ U dz y dz R dz T dz 

Assuming that R = constant and y = constant and multiplying the 
equation by dz result in 

dM _ dU 1_ dT 

M U " 2 T 


and rearranging yields 

dU dM 1 dT 
U _ M + 2 T 


Substituting the previously derived expression for dT/T yields 


52 = M + i (1 . r ,M» 52 

U M 2 ' 7 U 


dU r, 1 ,, > --.I dM 

— [1 --(1 - v ) . 


dU 


U ^1 - M 2 ^M 


dM 


The first term in equation (121) finally results in 



(1 + 6*/0) dU 1 + 6% 

dM 

U dz - ^ ) 

-- d Z 

M 


The second term of equation (121) can be expressed by 

_1_ d(pu) = 1 dU i dp 
pU dz U dz + p dz 
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The first right-hand term has already been expressed by the Maeh 
number expression. The second right-hand term is derived from the ideal 
gas equation as follows. 



Differentiating by parts in logarithmic form yields 


d(ln p ) - 

d (In P) 

- <1 (111 R) 

- d (In T) 

d(ln p) _ d(ln P) 

d(ln R) 

d(ln T) 

dz 

(1 7, 

dz 

" dz 


d(ln p) dp d(ln P) dP 

d (In R ) dR 

d(ln T) dT 

dp dz dP dz 

dR dz 

dT dz 

J. dp _ 1_ dP 1 dll 

1_ dT 


p dz P dz R dz 

T dz 



with R = constant and multiplication by dz results in 

d£ dP dT 

P P ' T 

Another expression for dp/p can be derived from the thermodynamic 
adiabatic relationship 

— — = const 


Differentiating by parts in logarithmic form yields 
d(ln P) - y d(ln p) d(ln const) 


Assuming that the specific heat ratio y — constant and dividing by 


dz result in 



d(ln P) 

- v 

d(ln p) d(ln const) 

dz 

i 

dz dz ’ 

d (hi P) 

dP 

d(ln p) dp 

dP 

dz 

^ dp dz ’ 
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0 


J_ op 1 dp 

P dz • ^ p dz 


1 



Final reari’angement results in 



Substitution yields 

dp _ dp dT 
p ~ Y p “ T 


or 


dp 

P 


(1 - y) = - 


dT 

T 


The previously derived temperature-Mach number relationship yields 


dp 

P 


Thus 


(1 - y) M 2 ^ 


dp = M 2 
P U 


Also derived was the velocity-Mach number relationship: 



1 

1 - 2 ~ -^ j M 2 M 


dM 


Therefore the second term in equation (121) yields 
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